diff --git a/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_Rev26.docx b/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_Rev26.docx new file mode 100644 index 0000000000..291484f998 Binary files /dev/null and b/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_Rev26.docx differ diff --git a/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_WakeExtentAndBuffering_Rev3.docx b/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_WakeExtentAndBuffering_Rev3.docx new file mode 100644 index 0000000000..5bc7691131 Binary files /dev/null and b/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_WakeExtentAndBuffering_Rev3.docx differ diff --git a/docs/OtherSupporting/FAST.Farm_Plan_Rev25.doc b/docs/OtherSupporting/FAST.Farm_Plan_Rev25.doc deleted file mode 100644 index 58ad2aad71..0000000000 Binary files a/docs/OtherSupporting/FAST.Farm_Plan_Rev25.doc and /dev/null differ diff --git a/docs/source/user/fast.farm/InputFiles.rst b/docs/source/user/fast.farm/InputFiles.rst index bcc524f4e7..223410552a 100644 --- a/docs/source/user/fast.farm/InputFiles.rst +++ b/docs/source/user/fast.farm/InputFiles.rst @@ -827,13 +827,15 @@ but include leading zeros. **WrDisDT** [sec] specifies the time step (inverse of the frame rate) of all disturbed wind data output files and must be an integer multiple -larger than or equal to **DT_Low**. This input is unused when -**WrDisWind** = FALSE and when **NOutDisWindXY**, **NOutDisWindYZ**, and -**NOutDisWindXZ** are set to zero. If the DEFAULT keyword is specified -in place of a numerical value, **WrDisDT** is set to **DT_Low**. Note -that the full high-resolution disturbed wind data output files are not -output at a frame rate of 1/**DT_High**, but are only output every -**WrDisDT** seconds. +larger than or equal to **DT_Low**. **WrDisDT** also sets the output +time step of the ambient wind and array effects wake-plane files written +to ``vtk_ff/wakes`` when **OutAllPlanes** = TRUE. This input is unused when +**WrDisWind** = FALSE, **OutAllPlanes** = FALSE, and **NOutDisWindXY**, +**NOutDisWindYZ**, and **NOutDisWindXZ** are set to zero. If the DEFAULT +keyword is specified in place of a numerical value, **WrDisDT** is set to +**DT_Low**. Note that the full high-resolution disturbed wind data output +files are not output at a frame rate of 1/**DT_High**, but are only output +every **WrDisDT** seconds. Visualizing the ambient wind and wake interactions can be useful for interpreting results and debugging problems. However, FAST.Farm will @@ -886,6 +888,10 @@ details on time-series results files. **OutAllPlanes** [-] Output all wake planes in VTK at all time steps. +This controls both the Wake Dynamics plane files written to ``vtk_ff_planes`` and the +wake-plane files written by the ambient wind and array effects module to ``vtk_ff/wakes`` +(see :numref:`FF:Output:Planes`). The ``vtk_ff/wakes`` files are written every **WrDisDT** +seconds, so increasing **WrDisDT** reduces the number of files written. Note: this option requires intensive writing to disk and will drastically slow down the simulation. DEFAULT is False. @@ -1134,24 +1140,49 @@ The sub-volume directories must follow the naming convention and ``1`` through **NumTurbines** for the high-resolution domain of each wind turbine. -- ** is a zero-padded integer index with the same field width as - **DirStartIndex**. The index for time step :math:`n` is - :math:`\text{DirStartNum} + n \times \Delta_\text{index}`, where - :math:`\text{DirStartNum}` is the integer value of **DirStartIndex** - and :math:`\Delta_\text{index}` is the stride between successive - directory indices (determined automatically by FAST.Farm by scanning - the available directories). - -During initialization, FAST.Farm reads the header of the starting -sub-volume for each domain (low-resolution sub-volume 0 and each -high-resolution sub-volume 1 through **NumTurbines**) to obtain the -grid properties, then calls the directory-discovery routine to confirm -that a sufficient number of time steps exist and that the grid -properties are consistent across all time steps. Specifically, FAST.Farm -verifies: - -- That at least **NumDT** low-resolution and high-resolution - directories are available for the simulation duration. +- ** is an integer index, zero-padded to at least the field width + of **DirStartIndex** and widened as needed when the step counter grows + past that width. + +FAST.Farm does **not** assume a constant stride between successive +directory indices. During initialization it reads the ``Header`` of each +directory matching the prefix and assigns it to a time step by the +simulation time recorded there: the directory used for time step +:math:`n` is the one whose header time is +:math:`t_\text{start} + n\,\Delta t`, where :math:`t_\text{start}` is the +header time of the directory named by **DirStartIndex** and +:math:`\Delta t` is **DT_Low-AMReX** for the low-resolution domain or +**DT_High-AMReX** for the high-resolution domains. + +This supports precursor data written with a varying solver time step -- +for example an AMR-Wind run that transitions from ``time.initial_dt`` to +``fixed_dt``, where the stride between directory indices changes but the +output interval in time does not -- provided the output interval itself +is uniform. + +Times are matched to within a tolerance of +:math:`\max(10^{-6}\ \mathrm{s},\ 10^{-9}(|t_\text{start}| + |t|))`, which +accommodates the floating-point drift accumulated by a long precursor +simulation while remaining far tighter than a single output interval. + +Having obtained the grid properties from the starting sub-volume of each +domain (low-resolution sub-volume 0 and each high-resolution sub-volume 1 +through **NumTurbines**), FAST.Farm verifies: + +- That every required time step is represented by a directory: steps 0 + through **NumDT** - 1 for the low-resolution domain, and steps 0 + through (**NumDT** - 1) x (**DT_Low-AMReX** / **DT_High-AMReX**) for each + high-resolution domain. A time step matched by no directory (missing + data) is a fatal error naming the time step and its expected simulation + time, as is a time step matched by more than one of the directories + scanned (for example, overlapping output left behind by a restart). The + duplicate check is not exhaustive: the scan stops early once every + required step has been matched, so a stale duplicate at a higher + directory index than the first directory past the simulation window is + not read. Directories whose header time falls between two required + steps are ignored, so the LES may write its sub-volumes at a finer + cadence than FAST.Farm reads them, provided the FAST.Farm time step is + an integer multiple of the output interval. - That the grid dimensions, origin, and spacing are identical across all time steps for a given sub-volume. @@ -1159,6 +1190,12 @@ verifies: - That the grid dimensions are the same for all high-resolution sub-volumes (one per turbine). +- That every high-resolution sub-volume is written at the same + simulation times as sub-volume 1. FAST.Farm uses a single + **DT_High-AMReX** for the whole farm, so a sub-volume on a different + set of times would supply its turbine a different instant than the + rest of the farm. + Because the grid properties are determined from the precursor data files, the user does not specify low- or high-resolution grid dimensions or origins in the FAST.Farm primary input file when using diff --git a/docs/source/user/fast.farm/OutputFiles.rst b/docs/source/user/fast.farm/OutputFiles.rst index 4e65e3accf..c63686b4f2 100644 --- a/docs/source/user/fast.farm/OutputFiles.rst +++ b/docs/source/user/fast.farm/OutputFiles.rst @@ -98,6 +98,8 @@ Wake dynamics Plane Files Setting the option **OutAllPlanes** to true in the main FAST.Farm input file will result in the wake planes of the Wake Dynamics module to be written. +The same option also controls the wake-plane VTK output of the ambient wind and array +effects module, written to the ``wakes`` subfolder of ``vtk_ff``. This option requires intensive writing to disk and will drastically slow down the simulation. The wake planes are written in VTK format, in the folder `vtk_ff_planes` at the root of the simulation directory. diff --git a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 index 04d44ae139..1d39ef8102 100644 --- a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 +++ b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 @@ -222,7 +222,11 @@ SUBROUTINE Farm_Initialize( farm, InputFile, ErrStat, ErrMsg ) call AllocAry( farm%p%MaxNumPlanes, farm%p%NumTurbines, 'farm%p%MaxNumPlanes', ErrStat2, ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName); if (Failed()) return do i=1,farm%p%NumTurbines ! Eventually, we will have different settings for different rotors - farm%p%MaxNumPlanes(i) = ceiling( 15.0 * ( WD_InitInput%InputFileData%NumDFull + WD_InitInput%InputFileData%NumDBuff ) / AWAE_InitInput%InputFileData%C_Meander ) + if (WD_InitInput%InputFileData%Mod_Wake == Mod_Wake_Polar) then + farm%p%MaxNumPlanes(i) = ceiling( 18.0 * ( WD_InitInput%InputFileData%NumDFull + WD_InitInput%InputFileData%NumDBuff ) / AWAE_InitInput%InputFileData%C_Meander ) + else + farm%p%MaxNumPlanes(i) = ceiling( 54.0 * ( WD_InitInput%InputFileData%NumDFull + WD_InitInput%InputFileData%NumDBuff ) ) + endif farm%p%MaxNumPlanes(i) = max( 2, min( farm%p%MaxNumPlanes(i) , farm%p%n_TMax + 2 ) ) end do @@ -245,6 +249,7 @@ SUBROUTINE Farm_Initialize( farm, InputFile, ErrStat, ErrMsg ) AWAE_InitInput%InputFileData%dt_low = farm%p%dt_low AWAE_InitInput%InputFileData%NumTurbines = farm%p%NumTurbines AWAE_InitInput%InputFileData%NumRadii = WD_InitInput%InputFileData%NumRadii + AWAE_InitInput%InputFileData%OutAllPlanes = WD_InitInput%InputFileData%OutAllPlanes AWAE_InitInput%MaxPlanes = MAXVAL(farm%p%MaxNumPlanes) AWAE_InitInput%InputFileData%WindFilePath = farm%p%WindFilePath AWAE_InitInput%n_high_low = farm%p%n_high_low @@ -698,6 +703,12 @@ SUBROUTINE Farm_InitWD( farm, WD_InitInp, ErrStat, ErrMsg ) WD_InitInp%TurbNum = nt WD_InitInp%MaxNumPlanes = farm%p%MaxNumPlanes(nt) WD_InitInp%OutFileRoot = farm%p%OutFileRoot + WD_InitInp%LowResBounds(1,1) = farm%p%X0_low + WD_InitInp%LowResBounds(2,1) = farm%p%Y0_low + WD_InitInp%LowResBounds(3,1) = farm%p%Z0_low + WD_InitInp%LowResBounds(1,2) = farm%p%X0_low + farm%p%dX_low * real(farm%p%nX_low - 1, ReKi) + WD_InitInp%LowResBounds(2,2) = farm%p%Y0_low + farm%p%dY_low * real(farm%p%nY_low - 1, ReKi) + WD_InitInp%LowResBounds(3,2) = farm%p%Z0_low + farm%p%dZ_low * real(farm%p%nZ_low - 1, ReKi) ! note that WD_Init has Interval as INTENT(IN) so, we don't need to worry about overwriting farm%p%dt_low here: call WD_Init( WD_InitInp, farm%WD(nt)%u, farm%WD(nt)%p, farm%WD(nt)%x, farm%WD(nt)%xd, farm%WD(nt)%z, & diff --git a/modules/awae/CMakeLists.txt b/modules/awae/CMakeLists.txt index 2079ef67a0..c7f42a576f 100644 --- a/modules/awae/CMakeLists.txt +++ b/modules/awae/CMakeLists.txt @@ -26,6 +26,7 @@ add_library(awaelib STATIC src/AWAE.f90 src/AWAE_IO.f90 src/AWAE_Types.f90 + src/AWAE_vtk.f90 src/amrex_utils.F90 ) target_link_libraries(awaelib awaelib_c ifwlib nwtclibs) diff --git a/modules/awae/src/AWAE.f90 b/modules/awae/src/AWAE.f90 index 8b329cd772..387952c567 100644 --- a/modules/awae/src/AWAE.f90 +++ b/modules/awae/src/AWAE.f90 @@ -27,6 +27,7 @@ module AWAE use NWTC_Library use AWAE_Types use AWAE_IO + use AWAE_vtk use InflowWind use IfW_FlowField use KdTree @@ -73,6 +74,12 @@ module AWAE contains +!---------------------------------------------------------------------------------------------------------------------------------- +!> Extract a 2D slice from a 3D vector field `V` at the requested physical coordinate `s` along the slice-normal axis. +!! The slice orientation is selected by `sliceType` (XYSlice, YZSlice, or XZSlice). The location `s` (in meters) is +!! converted to grid units using the grid origin `s0` and spacing `ds`, the two bracketing grid planes are identified, +!! and the output `slice` is filled by linear interpolation between them. If `s` coincides with the last grid index, +!! the upper bracketing index is clamped so no out-of-bounds access occurs. subroutine ExtractSlice( sliceType, s, s0, szs, sz1, sz2, ds, V, slice) integer(IntKi), intent(in ) :: sliceType !< Type of slice: XYSlice, YZSlice, XZSlice @@ -120,8 +127,16 @@ subroutine ExtractSlice( sliceType, s, s0, szs, sz1, sz2, ds, V, slice) end subroutine ExtractSlice !---------------------------------------------------------------------------------------------------------------------------------- -!> This subroutine -!! +!> Precompute, for every pair of adjacent wake planes (np, np+1) of every turbine, the geometric quantities that +!! describe the relative orientation of the two planes. For each pair, this routine evaluates the cosine and sine of +!! the angle between the plane normals `u%xhat_plane(:,np,nt)` and `u%xhat_plane(:,np+1,nt)` and uses them, together +!! with the offset between the plane centers `u%p_plane`, to determine whether the planes are (numerically) parallel. +!! When they are not parallel, the routine computes and caches in the misc-var struct `m` the perpendicular distances +!! from each plane center to the line of intersection of the two planes (`r_s`, `r_e`), the in-plane unit vectors +!! pointing from that intersection line toward each plane center (`rhat_s`, `rhat_e`), and the closest points on the +!! intersection line to each plane center (`pvec_cs`, `pvec_ce`). The boolean `m%parallelFlag(np,nt)` records the +!! parallel/non-parallel decision. These cached quantities are reused downstream (e.g., in `interp_planes_2_point`) to +!! interpolate the wake-plane center and orientation between adjacent skewed wake planes. subroutine ComputeLocals(n, u, p, y, m, errStat, errMsg) integer(IntKi), intent(in ) :: n !< Current simulation time increment (zero-based) type(AWAE_InputType), intent(in ) :: u !< Inputs at Time t @@ -830,8 +845,8 @@ subroutine LowResGridCalcOutput(n, u, p, xd, y, m, errStat, errMsg) ! - no messages if inside bounds, so put error handling inside if - call PlaneOutOfDomain(u%D_wake(np,nt),u%p_plane(:,np,nt),y%V_plane(:,np,nt),m%planeDomainExit(np,nt),ErrStat2,ErrMsg2) - if (m%planeDomainExit(np,nt) /= 0_IntKi) then + call PlaneOutOfDomain(p%y(p%NumRadii-1)+p%y(1),u%p_plane(:,np,nt),y%V_plane(:,np,nt),m%planeDomainExit(:,np,nt),ErrStat2,ErrMsg2) + if (any(m%planeDomainExit(:,np,nt) /= 0_IntKi)) then call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) cycle endif @@ -891,97 +906,72 @@ end function Failed !> Check if the center of this wwake plane has left the domain. !! If a plane exits the domain, or previously exited the domain: !! - Set warning about first time this plane leaves. - !! - Set component perpendicular to plane exit direction to kick it outside the domain entirely + !! - For each dimension independently, set velocity component to kick it outside the domain !! - Target distance outside boundary = D. Use a quadratic asymptotic distance per step to approach target distance. - !! - Add background flow in X or Y to keep the plane moving with others parallel to boundary it crossed (only using X and Y velocity) - !! NOTE: using m%planeDomainExit to track which boundary a plane crossed. - !! 0: Still in domain - !! +/-1: +/-X - !! +/-2: +/-Y - !! +/-3: +/-Z + !! - Add background flow to keep the plane drifting + !! NOTE: using m%planeDomainExit(3) to track which boundary a plane crossed in each dimension. + !! planeDomainExit(i): 0 = still in domain, -1 = crossed lower bound, +1 = crossed upper bound !! To understand intent, consider 2 cases for mean velocity in +X direction: !! plane exits +Y boundary: - !! 1. plane with get a kick towards one wake diameter outside +Y boundary + !! 1. plane will get a kick towards one wake diameter outside +Y boundary !! 2. overall farm velocity added to keep plane drifting in +X following the target Y location (some jitter due to farm level Y velocity term) !! plane exits +X boundary (travels beyond domain end in direction of overall flow) !! 1. plane will get a kick outside the end of the domain towards +X boundary plus wake diameter !! 2. farm velocity added will keep trying to push this plane further downstream, but step 1. will try to force it back. - !! --> effectively 1. and 2. will constant be working against each other to hold the plane somewhere near the target location beyond +X boundary, + !! --> effectively 1. and 2. will constantly be working against each other to hold the plane somewhere near the target location beyond +X boundary, !! but this shouldn't really matter as the plane will get dropped at some point. Even if multiple planes end up there, it shouldn't affect !! any planes still in bounds -- so we really don't care if it jitters around at all - subroutine PlaneOutOfDomain(D_Wake,p_plane,V_plane,planeDomainExit,ErrStat3,ErrMsg3) - real(ReKi), intent(in ) :: D_wake !< u%D_wake(np,nt) + subroutine PlaneOutOfDomain(PlaneHalfWidth,p_plane,V_plane,planeDomainExit,ErrStat3,ErrMsg3) + real(ReKi), intent(in ) :: PlaneHalfWidth !< half-width of wake plane + dr real(ReKi), intent(in ) :: p_plane(3) !< u%p_plane(:,np,nt) real(ReKi), intent(inout) :: V_plane(3) !< y%V_plane(:,np,nt) - integer(IntKi), intent(inout) :: planeDomainExit !< m%planeDomainExit(np,nt) + integer(IntKi), intent(inout) :: planeDomainExit(3) !< m%planeDomainExit(:,np,nt) per-dimension flag integer(IntKi), intent( out) :: ErrStat3 !< Error status of the operation character(ErrMsgLen), intent( out) :: ErrMsg3 !< Error message if errStat /= ErrID_None character(12) :: tmpStr12 !< for constructing error message real(ReKi) :: D_tgt !< target distance outside bounds - ! Step 1: did a plane that was in the low res domain just cross out? - ! If plane crossed boundary, set message and tracking of it - if (planeDomainExit == 0_IntKi) then - if (p_plane(1) < p%LowRes%oXYZ(1)) then ! lower x boundary - ErrStat3 = ErrID_Warn - tmpStr12 = 'lower-most X' - planeDomainExit = -1 - elseif ( p_plane(1) > p%LowRes%oXYZ(1) + p%LowRes%Size(1)) then ! upper x boundary - ErrStat3 = ErrID_Warn - tmpStr12 = 'upper-most X' - planeDomainExit = 1 - elseif ( p_plane(2) < p%LowRes%oXYZ(2)) then ! lower y boundary - ErrStat3 = ErrID_Warn - tmpStr12 = 'lower-most Y' - planeDomainExit = -2 - elseif ( p_plane(2) > p%LowRes%oXYZ(2) + p%LowRes%Size(2)) then ! upper y boundary - ErrStat3 = ErrID_Warn - tmpStr12 = 'upper-most Y' - planeDomainExit = 2 - elseif ( p_plane(3) < p%LowRes%oXYZ(3)) then ! lower z boundary - ErrStat3 = ErrID_Warn - tmpStr12 = 'lower-most Z' - planeDomainExit = -3 - elseif ( p_plane(3) > p%LowRes%oXYZ(3) + p%LowRes%Size(3)) then ! upper z boundary - ErrStat3 = ErrID_Warn - tmpStr12 = 'upper-most Z' - planeDomainExit = 3 - endif - if (errStat3 == ErrID_Warn) then - ErrMsg3 = 'The center of wake plane #'//trim(num2lstr(np))//' for turbine #'//trim(num2lstr(nt))//' has passed the ' & - //tmpStr12//' boundary of the low-resolution domain. Further warnings are suppressed.' + integer(IntKi) :: iDim !< loop counter over dimensions + character(1), parameter :: dimLabels(3) = (/'X','Y','Z'/) + + ! Step 1: did a plane that was in the low res domain just cross out in any dimension? + ! If plane crossed boundary, set message and tracking of it (check each dimension independently) + do iDim = 1, 3 + if (planeDomainExit(iDim) == 0_IntKi) then + if (p_plane(iDim) < p%LowRes%oXYZ(iDim)) then + ErrStat3 = ErrID_Warn + tmpStr12 = 'lower-most '//dimLabels(iDim) + planeDomainExit(iDim) = -1_IntKi + elseif (p_plane(iDim) > p%LowRes%oXYZ(iDim) + p%LowRes%Size(iDim)) then + ErrStat3 = ErrID_Warn + tmpStr12 = 'upper-most '//dimLabels(iDim) + planeDomainExit(iDim) = 1_IntKi + endif endif + end do + if (errStat3 == ErrID_Warn) then + ErrMsg3 = 'The center of wake plane #'//trim(num2lstr(np))//' for turbine #'//trim(num2lstr(nt))//' has passed the ' & + //tmpStr12//' boundary of the low-resolution domain. Further warnings are suppressed.' endif - - ! Step 2: for planes outside boundary (including one that just crossed outside) set velocity component to approach target offset. - ! asymptotically approach a distance D_wake away from the boundary (quadratic approach) - ! example: V at -Y boundary: - ! Vy = (Y_target - Y_pos) / (2 * DT) - select case (planeDomainExit) - case (0_IntKi) - return - case (-1_IntKi) ! Crossed -X - D_tgt = p%LowRes%oXYZ(1) - D_wake - V_plane(1) = (D_tgt - p_plane(1)) / (2.0_ReKi * real(p%dt_low,ReKi)) ! push towards (-X_bound - D_wake) - case ( 1_IntKi) ! Crossed +X - D_tgt = p%LowRes%oXYZ(1) + p%LowRes%Size(1) + D_wake - V_plane(1) = (D_tgt - p_plane(1)) / (2.0_ReKi * real(p%dt_low,ReKi)) ! push towards (+X_bound + D_wake) - case (-2_IntKi) ! Crossed -Y - D_tgt = p%LowRes%oXYZ(2) - D_wake - V_plane(2) = (D_tgt - p_plane(2)) / (2.0_ReKi * real(p%dt_low,ReKi)) ! push towards (-Y_bound - D_wake) - case ( 2_IntKi) ! Crossed +Y - D_tgt = p%LowRes%oXYZ(2) + p%LowRes%Size(2) + D_wake - V_plane(2) = (D_tgt - p_plane(2)) / (2.0_ReKi * real(p%dt_low,ReKi)) ! push towards (-Y_bound - D_wake) - case (-3_IntKi) ! Crossed -Z - D_tgt = p%LowRes%oXYZ(3) - D_wake - V_plane(3) = (D_tgt - p_plane(3)) / (2.0_ReKi * real(p%dt_low,ReKi)) ! push towards (-Z_bound - D_wake) - case ( 3_IntKi) ! Crossed +Z - D_tgt = p%LowRes%oXYZ(3) + p%LowRes%Size(3) + D_wake - V_plane(3) = (D_tgt - p_plane(3)) / (2.0_ReKi * real(p%dt_low,ReKi)) ! push towards (+Z_bound + D_wake) - end select - ! Step 3: add background XYZ flow to keep plane drifting (will have already returned on any planes still in bounds) + ! If still fully in domain, nothing to do + if (all(planeDomainExit == 0_IntKi)) return + + ! Step 2: add background XYZ flow to keep plane drifting V_plane(1:3) = V_plane(1:3) + xd%Ufarm(1:3) + ! Step 3: for each dimension where the plane is outside the boundary, set velocity component to reach target in one timestep. + ! example: V at -Y boundary: + ! Vy = (Y_target - Y_pos) / DT + do iDim = 1, 3 + if (planeDomainExit(iDim) == -1_IntKi) then + D_tgt = p%LowRes%oXYZ(iDim) - PlaneHalfWidth + V_plane(iDim) = (D_tgt - p_plane(iDim)) / real(p%dt_low,ReKi) + elseif (planeDomainExit(iDim) == 1_IntKi) then + D_tgt = p%LowRes%oXYZ(iDim) + p%LowRes%Size(iDim) + PlaneHalfWidth + V_plane(iDim) = (D_tgt - p_plane(iDim)) / real(p%dt_low,ReKi) + endif + end do + end subroutine PlaneOutOfDomain end subroutine LowResGridCalcOutput @@ -1186,7 +1176,7 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO integer(IntKi), intent( out) :: errStat !< Error status of the operation character(*), intent( out) :: errMsg !< Error message if errStat /= ErrID_None - character(1024) :: rootDir, baseName, OutFileVTKDir ! Simulation root dir, basename for outputs + character(1024) :: rootDir, baseName, OutFileVTKDir, OutFileVTKwakeDir ! Simulation root dir, basename for outputs integer(IntKi) :: i,j,nt,c ! loop counter real(ReKi) :: gridLoc ! Location of requested output slice in grid coordinates [0,sz-1] integer(IntKi) :: errStat2 ! temporary error status of the operation @@ -1240,7 +1230,6 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO ! AMReX Wind Parameters p%DirStartIndex = InitInp%InputFileData%DirStartIndex p%DirIndexLen = len_trim(InitInp%InputFileData%DirStartIndex) - read(p%DirStartIndex, *) p%DirStartNum ! Wake Added Turbulence (WAT) Parameters p%WAT_Enabled = InitInp%WAT_Enabled @@ -1267,13 +1256,35 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO ! --- Vtk Outputs call GetPath( p%OutFileRoot, rootDir, baseName ) - OutFileVTKDir = trim(rootDir) // 'vtk_ff' ! Directory for VTK outputs - p%OutFileVTKRoot = trim(rootDir) // 'vtk_ff' // PathSep // trim(baseName) ! Basename for VTK files + OutFileVTKDir = trim(rootDir) // 'vtk_ff' ! Directory for VTK outputs + p%OutFileFFvtkRoot = trim(OutFileVTKDir) // PathSep // trim(baseName) ! Basename for VTK files p%VTK_tWidth = CEILING( log10( real(p%NumDT, ReKi)/real(p%WrDisSkp1, ReKi) ) + 1) ! Length for time stamp if (p%WrDisWind .or. p%NOutDisWindXY>0 .or. p%NOutDisWindYZ>0 .or. p%NOutDisWindXZ>0) then - call MKDIR(OutFileVTKDir) ! creating output directory + call MKDIR(OutFileVTKDir) end if + ! Wake-plane VTK output is controlled solely by the OutAllPlanes input flag: it is + ! IO-intensive, so it must not be implied by any other VTK output request. + p%WrPlanes = InitInp%InputFileData%OutAllPlanes + + ! Setup wake plane writing + if (p%WrPlanes) then + OutFileVTKwakeDir = trim(OutFileVTKDir) // PathSep // 'wakes' ! Directory for VTK wake outputs + call MKDIR(OutFileVTKDir) ! we may not be writing out any other vtk, so create dir if doesn't exist + call MKDIR(OutFileVTKwakeDir) + p%OutFileFFvtkWakeRoot = trim(OutFileVTKwakeDir) // PathSep // trim(baseName) ! Basename for VTK wake files + p%VTK_tWidthPlanes = CEILING( log10(real(max(p%MaxPlanes, 1), ReKi)) + 1) ! length for the number of planes + ! Since a huge number of wake planes will be empty initially, we don't want to write all those out. + p%OutFileFFvtkWakeNullData = "FF.WakePlane_Null.vtk" + + ! track when a plane is first written out + allocate( m%WakeVTK_StartN(0:p%MaxPlanes-1,p%NumTurbines), stat=ErrStat2); if (Failed0('Could not allocate memory for m%WakeVTK_StartN')) return; + m%WakeVTK_StartN = huge(1_IntKi) + + ! write out the null wake plane + call Write_NullPlane(OutFileVTKwakeDir, p) + endif + ! Plane grids allocate( p%y(-p%Numradii+1:p%NumRadii-1), stat=errStat2); if (Failed0('Could not allocate memory for p%y.')) return; allocate( p%z(-p%Numradii+1:p%NumRadii-1), stat=errStat2); if (Failed0('Could not allocate memory for p%z.')) return; @@ -1442,10 +1453,14 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO allocate ( u%D_wake ( 0:p%MaxPlanes-1,1:p%NumTurbines), STAT=ErrStat2 ); if (Failed0('u%D_wake.' )) return; allocate ( u%WAT_k (1-p%NumRadii:p%NumRadii-1, 1-p%NumRadii:p%NumRadii-1, 0:p%MaxPlanes-1,1:p%NumTurbines), STAT=ErrStat2 ); if (Failed0('u%WAT_k.' )) return; - u%NumPlanes = 2.0_ReKi - u%Vx_wake=0.0_ReKi - u%Vy_wake=0.0_ReKi - u%Vz_wake=0.0_ReKi + u%NumPlanes = 2.0_ReKi + u%xhat_plane = 0.0_ReKi + u%p_plane = 0.0_ReKi + u%Vx_wake = 0.0_ReKi + u%Vy_wake = 0.0_ReKi + u%Vz_wake = 0.0_ReKi + u%D_wake = 0.0_ReKi + u%WAT_k = 0.0_ReKi !---------------------------------------------------------------------------- @@ -1471,9 +1486,9 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO end do ! This next step is not strictly necessary - y%V_plane = 0.0_Reki - y%Vx_wind_disk = 0.0_Reki - y%TI_amb = 0.0_Reki + y%V_plane = 0.0_Reki + y%Vx_wind_disk = 0.0_Reki + y%TI_amb = 0.0_Reki !---------------------------------------------------------------------------- ! Initialize misc @@ -1485,12 +1500,15 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO if ( p%NOutDisWindXY > 0 ) then ALLOCATE ( m%OutVizXYPlane(3,p%LowRes%nXYZ(1), p%LowRes%nXYZ(2),1) , STAT=ErrStat2 ); if (Failed0('the Fast.Farm OutVizXYPlane arrays.')) return; + m%OutVizXYPlane = 0.0_SiKi end if if ( p%NOutDisWindYZ > 0 ) then ALLOCATE ( m%OutVizYZPlane(3,p%LowRes%nXYZ(2), p%LowRes%nXYZ(3),1) , STAT=ErrStat2 ); if (Failed0('the Fast.Farm OutVizYZPlane arrays.')) return; + m%OutVizYZPlane = 0.0_SiKi end if if ( p%NOutDisWindXZ > 0 ) then ALLOCATE ( m%OutVizXZPlane(3,p%LowRes%nXYZ(1), p%LowRes%nXYZ(3),1) , STAT=ErrStat2 ); if (Failed0('the Fast.Farm OutVizXZPlane arrays.')) return; + m%OutVizXZPlane = 0.0_SiKi end if ! miscvars to avoid the allocation per timestep @@ -1498,20 +1516,25 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO allocate(m%Vamb_low( 3, 0:p%LowRes%nXYZ(1)-1 , 0:p%LowRes%nXYZ(2)-1 , 0:p%LowRes%nXYZ(3)-1 ), STAT=errStat2); if (Failed0('m%Vamb_low.' )) return; allocate(m%Vdist_low( 3, 0:p%LowRes%nXYZ(1)-1 , 0:p%LowRes%nXYZ(2)-1 , 0:p%LowRes%nXYZ(3)-1 ), STAT=errStat2); if (Failed0('m%Vdist_low.' )) return; allocate(m%Vdist_low_full( 3, 0:p%LowRes%nXYZ(1)-1 , 0:p%LowRes%nXYZ(2)-1 , 0:p%LowRes%nXYZ(3)-1 ), STAT=errStat2); if (Failed0('m%Vdist_low_full')) return; + m%Vamb_lowpol = 0.0_ReKi + m%Vamb_low = 0.0_SiKi + m%Vdist_low = 0.0_SiKi + m%Vdist_low_full = 0.0_SiKi allocate(m%Vamb_high(1:p%NumTurbines), STAT=ErrStat2); if (Failed0('Could not allocate memory for m%Vamb_high.')) return; do nt = 1, p%NumTurbines allocate(m%Vamb_high(nt)%data(3,0:p%HighRes(nt)%nXYZ(1)-1, 0:p%HighRes(nt)%nXYZ(2)-1, 0:p%HighRes(nt)%nXYZ(3)-1, 0:p%n_high_low_p1), STAT=ErrStat2) if (Failed0('m%Vamb_high%data.')) return; + m%Vamb_high(nt)%data = 0.0_SiKi end do - allocate(m%parallelFlag( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%parallelFlag.')) return; - allocate(m%r_s( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%r_s.' )) return; - allocate(m%r_e( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%r_e.' )) return; - allocate(m%rhat_s( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%rhat_s.' )) return; - allocate(m%rhat_e( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%rhat_e.' )) return; - allocate(m%pvec_cs( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%pvec_cs.' )) return; - allocate(m%pvec_ce( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%pvec_ce.' )) return; + allocate(m%parallelFlag( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%parallelFlag.')) return; m%parallelFlag = .false. + allocate(m%r_s( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%r_s.' )) return; m%r_s = 0.0_ReKi + allocate(m%r_e( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%r_e.' )) return; m%r_e = 0.0_ReKi + allocate(m%rhat_s( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%rhat_s.' )) return; m%rhat_s = 0.0_ReKi + allocate(m%rhat_e( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%rhat_e.' )) return; m%rhat_e = 0.0_ReKi + allocate(m%pvec_cs( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%pvec_cs.' )) return; m%pvec_cs = 0.0_ReKi + allocate(m%pvec_ce( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%pvec_ce.' )) return; m%pvec_ce = 0.0_ReKi ! WAT - store array of disk average velocities for all turbines call AllocAry(m%V_amb_low_disk,3,p%NumTurbines,'m%V_amb_low_disk', ErrStat2, ErrMsg2); if(Failed()) return; @@ -1540,9 +1563,9 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO ! Initialize the KdTree with no active points call kdtree_build(m%KdT, m%AllPlanePoints(:,1:1), n_max=p%MaxPlanes*p%NumTurbines) - ! track if a plan has left the domain (all planes start in domain). - ! Value indicates edge number (+/-1: +/-X, +/-2: +/-Y, +/-3: +/-Z) the plane crossed - allocate(m%planeDomainExit(0:p%MaxPlanes-1,1:p%NumTurbines), STAT=ErrStat2); if (Failed0('m%planeDomainExit.')) return; + ! track if a plane has left the domain (all planes start in domain). + ! Per-dimension flag: 0 = still in domain, -1 = crossed lower bound, +1 = crossed upper bound + allocate(m%planeDomainExit(3,0:p%MaxPlanes-1,1:p%NumTurbines), STAT=ErrStat2); if (Failed0('m%planeDomainExit.')) return; m%planeDomainExit = 0_IntKi ! Read-in the ambient wind data for the initial calculate output @@ -1674,6 +1697,26 @@ subroutine AWAE_End( u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) errStat = ErrID_None errMsg = "" + ! Write .vtk.series files for disturbed-wind output slices + do nt = 1, p%NOutDisWindXY + if (.not. p%OutDisWindZvalid(nt)) cycle + call Write_DisWind_Series(p, "DisXY", nt) + end do + do nt = 1, p%NOutDisWindYZ + if (.not. p%OutDisWindXvalid(nt)) cycle + call Write_DisWind_Series(p, "DisYZ", nt) + end do + do nt = 1, p%NOutDisWindXZ + if (.not. p%OutDisWindYvalid(nt)) cycle + call Write_DisWind_Series(p, "DisXZ", nt) + end do + + if (p%WrPlanes) then + ! Write final ParaView .vtk.series files for wake planes + call Write_WakePlane_Series(p, m) + call Write_WireFrame_Series(p) + endif + ! Destroy InflowWind data select case(p%Mod_AmbWind) case (2) @@ -2085,7 +2128,7 @@ subroutine AWAE_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg call ExtractSlice(XYSlice, p%OutDisWindZ(k), p%LowRes%oXYZ(3), p%LowRes%nXYZ(3), p%LowRes%nXYZ(1), p%LowRes%nXYZ(2), p%LowRes%dXYZ(3), m%Vdist_low_full, m%outVizXYPlane(:,:,:,1)) ! Create the output vtk file with naming /Low/DisXY.t.vtk - FileName = trim(p%OutFileVTKRoot)//".Low.DisXY"//PlaneNumStr//"."//trim(Tstr)//".vtk" + FileName = trim(p%OutFileFFvtkRoot)//".Low.DisXY"//PlaneNumStr//"."//trim(Tstr)//".vtk" call WrVTK_SP_header(FileName, "Low resolution, disturbed wind of XY Slice at time = "//trim(num2lstr(t))//" seconds.", Un, ErrStat2, ErrMsg2 ); if (Failed()) return; call WrVTK_SP_vectors3D(Un, "Velocity", & [p%LowRes%nXYZ(1), p%LowRes%nXYZ(2), 1_IntKi], & @@ -2102,7 +2145,7 @@ subroutine AWAE_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg call ExtractSlice(YZSlice, p%OutDisWindX(k), p%LowRes%oXYZ(1), p%LowRes%nXYZ(1), p%LowRes%nXYZ(2), p%LowRes%nXYZ(3), p%LowRes%dXYZ(1), m%Vdist_low_full, m%outVizYZPlane(:,:,:,1)) ! Create the output vtk file with naming /Low/DisYZ.t.vtk - FileName = trim(p%OutFileVTKRoot)//".Low.DisYZ"//PlaneNumStr//"."//trim(Tstr)//".vtk" + FileName = trim(p%OutFileFFvtkRoot)//".Low.DisYZ"//PlaneNumStr//"."//trim(Tstr)//".vtk" call WrVTK_SP_header(FileName, "Low resolution, disturbed wind of YZ Slice at time = "//trim(num2lstr(t))//" seconds.", Un, ErrStat2, ErrMsg2 ); if (Failed()) return; call WrVTK_SP_vectors3D(Un, "Velocity", & [1, p%LowRes%nXYZ(2), p%LowRes%nXYZ(3)], & @@ -2119,7 +2162,7 @@ subroutine AWAE_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg call ExtractSlice(XZSlice, p%OutDisWindY(k), p%LowRes%oXYZ(2), p%LowRes%nXYZ(2), p%LowRes%nXYZ(1), p%LowRes%nXYZ(3), p%LowRes%dXYZ(2), m%Vdist_low_full, m%outVizXZPlane(:,:,:,1)) ! Create the output vtk file with naming /Low/DisXZ.t.vtk - FileName = trim(p%OutFileVTKRoot)//".Low.DisXZ"//PlaneNumStr//"."//trim(Tstr)//".vtk" + FileName = trim(p%OutFileFFvtkRoot)//".Low.DisXZ"//PlaneNumStr//"."//trim(Tstr)//".vtk" call WrVTK_SP_header(FileName, "Low resolution, disturbed wind of XZ Slice at time = "//trim(num2lstr(t))//" seconds.", Un, ErrStat2, ErrMsg2); if (Failed()) return; call WrVTK_SP_vectors3D(Un, "Velocity", & [p%LowRes%nXYZ(1), 1, p%LowRes%nXYZ(3)], & @@ -2128,6 +2171,21 @@ subroutine AWAE_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg if (Failed()) return end do + if (p%WrPlanes) then + !------------------------------------------------------------------------- + ! Write a VTK polydata file containing the four corners of every active + ! wake plane for every turbine as a set of quads. + !------------------------------------------------------------------------- + call Write_Planes_WireFrame(p, u, t, Tstr) + + !------------------------------------------------------------------------- + ! Write one VTK STRUCTURED_GRID file per wake plane (per turbine) with + ! the wake velocity sampled on the plane's structured Y-Z grid, plus a + ! ParaView .vtk.series JSON index of all wake-plane files written so far. + !------------------------------------------------------------------------- + call Write_Planes_Data(p, u, m, n, t, Tstr) + endif + #ifdef FF_TIMING_PRINTS call AWAE_AddStageTiming('WriteDisWind', tmSer0, tmPar0) #endif @@ -2135,6 +2193,7 @@ subroutine AWAE_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg contains + logical function Failed() call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) Failed = ErrStat >= AbortErrLev @@ -2666,5 +2725,4 @@ real(ReKi) function testf(y,z) end function end subroutine - end module AWAE diff --git a/modules/awae/src/AWAE_IO.f90 b/modules/awae/src/AWAE_IO.f90 index 663c9e2dc5..1f6987cdb6 100644 --- a/modules/awae/src/AWAE_IO.f90 +++ b/modules/awae/src/AWAE_IO.f90 @@ -80,7 +80,7 @@ subroutine WriteDisWindFiles( n, WrDisSkp1, p, y, m, errStat, errMsg ) ! TimeStamp write(Tstr, '(i' // trim(Num2LStr(p%VTK_tWidth)) //'.'// trim(Num2LStr(p%VTK_tWidth)) // ')') n_out ! TODO use n instead.. - FileName = trim(p%OutFileVTKRoot)//".Low.Dis."//trim(Tstr)//".vtk" + FileName = trim(p%OutFileFFvtkRoot)//".Low.Dis."//trim(Tstr)//".vtk" call WrVTK_SP_header( FileName, "Low resolution disturbed wind for time = "//trim(num2lstr(t_out))//" seconds.", Un, errStat2, errMsg2 ) call SetErrStat(errStat2, errMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return @@ -92,7 +92,7 @@ subroutine WriteDisWindFiles( n, WrDisSkp1, p, y, m, errStat, errMsg ) ! We are only writing out the first of the high res data for a given low res time step ! NOTE: y%Vdist_high(nt)%data(:,:,:,:,1) is at T=t_low, and index 0 is at T=t_low-DT_high - FileName = trim(p%OutFileVTKRoot)//".HighT"//trim(num2lstr(nt))//".Dis."//trim(Tstr)//".vtk" + FileName = trim(p%OutFileFFvtkRoot)//".HighT"//trim(num2lstr(nt))//".Dis."//trim(Tstr)//".vtk" call WrVTK_SP_header( FileName, "High resolution disturbed wind for time = "//trim(num2lstr(t_out))//" seconds.", Un, errStat2, errMsg2 ) call SetErrStat(errStat2, errMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return @@ -171,24 +171,59 @@ subroutine ReadWindAMReX(sv, n, p, Vamb, ErrStat, ErrMsg) integer(IntKi), intent( out) :: ErrStat !< Error status of the operation character(*), intent( out) :: ErrMsg !< Error message if errStat /= ErrID_None + character(*), parameter :: RoutineName = 'ReadWindAMReX' character(len=2048) :: FileName ! Name of output file - character(len=12) :: DirIndex ! Directory index suffix - integer(IntKi) :: i + character(len=16) :: DirIndex ! Directory index suffix + character(len=16) :: NumStr ! Directory index without padding + integer(IntKi) :: DirIndexNum ! Directory index for this time step + integer(IntKi) :: nLo, nHi ! Bounds of the index table - ! If sub-volume is 0 then this is low-resolution file + ErrStat = ErrID_None + ErrMsg = "" + + ! Look up the directory index for this time step. Sub-volume 0 is the low-resolution domain, + ! sub-volume 1+ is turbine sv's high-resolution domain. Both tables are indexed by the 0-based + ! time step number and were built in AWAE_IO_InitGridInfo by matching header times, so they + ! carry no assumption that the LES advanced at a constant time step. if (sv == 0) then - write(DirIndex,'(i'//trim(Num2LStr(p%DirIndexLen))//')') p%DirStartNum + p%DirIndexDeltaLow * n + if (.not. allocated(p%DirIndexLow)) then + call SetErrStat(ErrID_Fatal, 'the low-resolution AMReX directory table was never populated; '// & + 'AWAE_IO_InitGridInfo did not complete successfully.', ErrStat, ErrMsg, RoutineName) + return + end if + nLo = lbound(p%DirIndexLow, 1); nHi = ubound(p%DirIndexLow, 1) + if (n < nLo .or. n > nHi) then + call SetErrStat(ErrID_Fatal, 'requested low-resolution AMReX time step '//trim(Num2LStr(n))// & + ', but only steps '//trim(Num2LStr(nLo))//' through '//trim(Num2LStr(nHi))// & + ' were located during initialization.', ErrStat, ErrMsg, RoutineName) + return + end if + DirIndexNum = p%DirIndexLow(n) else - write(DirIndex,'(i'//trim(Num2LStr(p%DirIndexLen))//')') p%DirStartNum + p%DirIndexDeltaHigh * n + if (.not. allocated(p%DirIndexHigh)) then + call SetErrStat(ErrID_Fatal, 'the high-resolution AMReX directory table was never populated; '// & + 'AWAE_IO_InitGridInfo did not complete successfully.', ErrStat, ErrMsg, RoutineName) + return + end if + nLo = lbound(p%DirIndexHigh, 1); nHi = ubound(p%DirIndexHigh, 1) + if (n < nLo .or. n > nHi) then + call SetErrStat(ErrID_Fatal, 'requested high-resolution AMReX time step '//trim(Num2LStr(n))// & + ' for sub-volume '//trim(Num2LStr(sv))//', but only steps '//trim(Num2LStr(nLo))// & + ' through '//trim(Num2LStr(nHi))//' were located during initialization.', & + ErrStat, ErrMsg, RoutineName) + return + end if + DirIndexNum = p%DirIndexHigh(n) end if - ! Prepend zeros in front of index number - do i = 1, p%DirIndexLen - if (DirIndex(i:i) /= " ") exit - DirIndex(i:i) = '0' - end do + ! Left-pad the index with zeros to DirIndexLen characters, widening the field when the + ! index needs more digits. DirIndexLen is a *minimum* width, matching amrex::Concatenate, + ! which pads with setw(mindigits); a fixed-width edit descriptor would overflow to '****' + ! once a run passes 10**DirIndexLen steps. + write(NumStr,'(I0)') DirIndexNum + DirIndex = repeat('0', max(0, p%DirIndexLen - len_trim(NumStr)))//trim(NumStr) - FileName = trim(p%WindFilePath)//"_"//trim(num2lstr(sv))//"_"//DirIndex(1:p%DirIndexLen) + FileName = trim(p%WindFilePath)//"_"//trim(num2lstr(sv))//"_"//trim(DirIndex) call amrex_read_data(FileName, Vamb, ErrStat, ErrMsg) end subroutine @@ -249,8 +284,9 @@ subroutine AWAE_IO_InitGridInfo(InitInp, p, InitOut, errStat, errMsg) integer(IntKi) :: nChunksX, nChunksY integer(IntKi) :: nChunkPointsX, nChunkPointsY integer(IntKi), allocatable :: ChunkIndicesX(:,:), ChunkIndicesY(:,:) - integer(IntKi) :: StartIndexNum, IndexDelta - + integer(IntKi) :: NumStepHigh ! Number of high-res time slices required + integer(IntKi), allocatable :: DirIndexTmp(:) ! Per-sub-volume index table, for cross-checking + errStat = ErrID_None errMsg = "" @@ -296,9 +332,9 @@ subroutine AWAE_IO_InitGridInfo(InitInp, p, InitOut, errStat, errMsg) call amrex_read_header(FileName, Time, dims, gridSpacing, origin, ErrStat2, ErrMsg2) if (Failed()) return - ! Search directory for time slices of this sub-volume - call amrex_find_subvols(p%WindFilePath, 0, p%dt_low, p%NumDT, p%DirStartIndex, & - StartIndexNum, p%DirIndexDeltaLow, ErrStat2, ErrMsg2) + ! Build the time step -> directory index table for the low-resolution sub-volume + call amrex_find_subvols(p%WindFilePath, 0, p%dt_low, p%NumDT, trim(p%DirStartIndex), & + p%DirIndexLow, ErrStat2, ErrMsg2) if (Failed()) return end select @@ -328,7 +364,6 @@ subroutine AWAE_IO_InitGridInfo(InitInp, p, InitOut, errStat, errMsg) p%LowRes%nXYZ = dims p%LowRes%nPoints = product(dims) p%LowRes%Size = gridSpacing * real(dims - 1, ReKi) - p%LowRes%Center = origin + 0.5_ReKi * p%LowRes%Size ! Polar data p%dPol = (gridSpacing(1)+gridSpacing(2)+gridSpacing(3))/3.0_ReKi @@ -500,18 +535,37 @@ subroutine AWAE_IO_InitGridInfo(InitInp, p, InitOut, errStat, errMsg) call amrex_read_header(FileName, Time, dims, gridSpacing, origin, ErrStat2, ErrMsg2) if (Failed()) return - ! Search directory for time slices of this sub-volume - call amrex_find_subvols(p%WindFilePath, nt, p%dt_high, p%NumDT*p%n_high_low-1, p%DirStartIndex, & - StartIndexNum, IndexDelta, ErrStat2, ErrMsg2) + ! Number of high-resolution time slices actually read by AWAE_UpdateStates: + ! indices n*n_high_low + i_hl for n = 0..NumDT-1 and i_hl = 0..n_high_low, where + ! i_hl is forced to 0 on the final low-res step. The largest index reached is + ! therefore (NumDT-1)*n_high_low, so (NumDT-1)*n_high_low + 1 slices are needed. + NumStepHigh = (p%NumDT - 1)*p%n_high_low + 1 + + ! Build the time step -> directory index table for this sub-volume + call amrex_find_subvols(p%WindFilePath, nt, p%dt_high, NumStepHigh, trim(p%DirStartIndex), & + DirIndexTmp, ErrStat2, ErrMsg2) if (Failed()) return - ! If first turbine, save index delta, otherwise ensure that it is the same + ! All high-resolution sub-volumes must be written at the same simulation times: FAST.Farm + ! has a single DT_High for the whole farm, so a sub-volume on a different set of times + ! would silently supply the wrong instant to its turbine. Keep sub-volume 1's table and + ! verify the rest match it element by element -- comparing only the stride, as was done + ! previously, accepts a sequence uniformly offset from the others. if (nt == 1) then - p%DirIndexDeltaHigh = IndexDelta - else if (p%DirIndexDeltaHigh /= IndexDelta) then - call SetErrStat(ErrID_Fatal, "got different index delta for sub-volume "//trim(Num2LStr(nt))//" than for sub-volume 1", & - ErrStat, ErrMsg, RoutineName) - return + call move_alloc(DirIndexTmp, p%DirIndexHigh) + else + do n = 0, NumStepHigh - 1 + if (DirIndexTmp(n) /= p%DirIndexHigh(n)) then + call SetErrStat(ErrID_Fatal, "high-resolution AMReX sub-volume "//trim(Num2LStr(nt))// & + " is not written at the same simulation times as sub-volume 1. At high-resolution"// & + " time step "//trim(Num2LStr(n))//" sub-volume "//trim(Num2LStr(nt))// & + " uses directory index "//trim(Num2LStr(DirIndexTmp(n)))//" but sub-volume 1 uses "// & + trim(Num2LStr(p%DirIndexHigh(n)))//". All sub-volumes must be output at the same times.", & + ErrStat, ErrMsg, RoutineName) + return + end if + end do + deallocate(DirIndexTmp) end if end select diff --git a/modules/awae/src/AWAE_Registry.txt b/modules/awae/src/AWAE_Registry.txt index 7a52d7e940..0ebb17df3b 100644 --- a/modules/awae/src/AWAE_Registry.txt +++ b/modules/awae/src/AWAE_Registry.txt @@ -40,6 +40,7 @@ typedef ^ ^ ReKi OutDisWindX {:} typedef ^ ^ IntKi NOutDisWindXZ - - - "Number of XZ planes for output of disturbed wind data across the low-resolution domain to /Low/DisXZ..t.vtk [0 to 9]" - typedef ^ ^ ReKi OutDisWindY {:} - - "Y coordinates of XZ planes for output of disturbed wind data across the low-resolution domain [1 to NOutDisWindXZ]" meters typedef ^ ^ DbKi WrDisDT - - - "The time between vtk outputs [must be a multiple of the low resolution time step]" s +typedef ^ ^ LOGICAL OutAllPlanes - .false. - "Output all wake planes in VTK at all time steps" - typedef ^ ^ LOGICAL ChkWndFiles - - - "Check all the ambient wind files for data consistency (flag)" - typedef ^ ^ IntKi Mod_Meander - - - "Spatial filter model for wake meandering {1: uniform, 2: truncated jinc, 3: windowed jinc} [DEFAULT=2]" - typedef ^ ^ ReKi C_Meander - - - "Calibrated parameter for wake meandering [>=1.0] [DEFAULT=1.9]" - @@ -127,13 +128,13 @@ typedef ^ MiscVarType IntKi iPlaneTurbTurb {:}{:}{:} - - "Fir typedef ^ MiscVarType IntKi iPlaneTurbChunk {:}{:}{:} - - "First and Last plane index by source turbine and destination chunk index" - typedef ^ MiscVarType Logical LowResChunkHasWake {:} - - "Low-res gridFirst and Last plane index by source turbine and destination chunk index" - typedef ^ MiscVarType ReKi MaxWakePointSep - - - "Maximum separation between wake points" - -typedef ^ MiscVarType Logical parallelFlag {:}{:} - - "" - -typedef ^ MiscVarType ReKi r_s {:}{:} - - "" - -typedef ^ MiscVarType ReKi r_e {:}{:} - - "" - -typedef ^ MiscVarType ReKi rhat_s {:}{:}{:} - - "" - -typedef ^ MiscVarType ReKi rhat_e {:}{:}{:} - - "" - -typedef ^ MiscVarType ReKi pvec_cs {:}{:}{:} - - "" - -typedef ^ MiscVarType ReKi pvec_ce {:}{:}{:} - - "" - +typedef ^ MiscVarType Logical parallelFlag {:}{:} - - "Flag indicating whether adjacent wake planes np and np+1 are (numerically) parallel for each turbine; dims: (plane-pair index np, turbine index nt)" - +typedef ^ MiscVarType ReKi r_s {:}{:} - - "Perpendicular distance from the start wake-plane center p_plane(:,np,nt) to the intersection line of adjacent (non-parallel) wake planes np and np+1; dims: (plane-pair index np, turbine index nt)" m +typedef ^ MiscVarType ReKi r_e {:}{:} - - "Perpendicular distance from the end wake-plane center p_plane(:,np+1,nt) to the intersection line of adjacent (non-parallel) wake planes np and np+1; dims: (plane-pair index np, turbine index nt)" m +typedef ^ MiscVarType ReKi rhat_s {:}{:}{:} - - "Unit vector lying in the start wake plane (np) pointing from the plane-plane intersection line toward the start plane center; dims: (XYZ component, plane-pair index np, turbine index nt)" - +typedef ^ MiscVarType ReKi rhat_e {:}{:}{:} - - "Unit vector lying in the end wake plane (np+1) pointing from the plane-plane intersection line toward the end plane center; dims: (XYZ component, plane-pair index np, turbine index nt)" - +typedef ^ MiscVarType ReKi pvec_cs {:}{:}{:} - - "Closest point on the plane-plane intersection line to the start wake-plane center, p_plane(:,np,nt) - r_s*rhat_s; dims: (XYZ component, plane-pair index np, turbine index nt)" m +typedef ^ MiscVarType ReKi pvec_ce {:}{:}{:} - - "Closest point on the plane-plane intersection line to the end wake-plane center, p_plane(:,np+1,nt) - r_e*rhat_e; dims: (XYZ component, plane-pair index np, turbine index nt)" m # typedef ^ MiscVarType SiKi outVizXYPlane {:}{:}{:}{:} - - "An array holding the output data for a 2D visualization slice" - typedef ^ MiscVarType SiKi outVizYZPlane {:}{:}{:}{:} - - "An array holding the output data for a 2D visualization slice" - @@ -145,7 +146,10 @@ typedef ^ MiscVarType InflowWind_OutputType y_IfW_High {:} - - "InflowWi #wake added turbulence typedef ^ MiscVarType ReKi V_amb_low_disk {:}{:} - - "Rotor averaged ambiend wind speed for each wind turbine (3 x nWT)" m/s -typedef ^ MiscVarType IntKi planeDomainExit {:}{:} 0 - "Value indicates edge number (0: still in domain, +/-1: +/-X, +/-2: +/-Y, +/-3: +/-Z) the plane crossed" - +typedef ^ MiscVarType IntKi planeDomainExit {:}{:}{:} 0 - "Per-dimension flag (0: still in domain, -1: crossed lower bound, +1: crossed upper bound) for each plane [dim,plane,turbine]" - + +#visualization +typedef ^ MiscVarType IntKi WakeVTK_StartN {:}{:} - - "VTK output index (n/WrDisSkp1) when wake plane is first written. Indices [wakenum,turbnum]" - # Low-resolution grid chunk data typedef ^ LRGChunkType IntKi iChunk {3} - - "XYZ index of chunk" - @@ -166,7 +170,6 @@ typedef ^ LRGParamType IntKi nXYZ {3} - - "Number typedef ^ LRGParamType IntKi nPoints - - - "Number of spatial nodes" - typedef ^ LRGParamType ReKi GridPoints {:}{:} - - "XYZ components (global positions) of the spatial discretization of the grid" m typedef ^ LRGParamType ReKi Size {3} - - "XYZ size of the grid" m -typedef ^ LRGParamType ReKi Center {3} - - "XYZ coordinates of the grid center" m typedef ^ LRGParamType LRGChunkType WakeChunks {:} - - "Chunks for updating grid from wake" - @@ -207,9 +210,8 @@ typedef ^ ParameterType ReKi C_ScaleDiam - - - "Normalize typedef ^ ParameterType IntKi Mod_Projection - - - "Switch to select how the wake plane velocity is projected in AWAE {1: keep all components, 2: project against plane normal} or DEFAULT [DEFAULT=1: if Mod_Wake is 1 or 3, or DEFAULT=2: if Mod_Wake is 2]" typedef ^ ParameterType character(12) DirStartIndex - - - "Starting directory index suffix for AMReX wind" - typedef ^ ParameterType IntKi DirIndexLen - - - "Number of characters in directory index" -typedef ^ ParameterType IntKi DirStartNum - - - "Starting directory index number for AMReX wind" -typedef ^ ParameterType IntKi DirIndexDeltaLow - - - "Directory index delta for low-resolution AMReX wind" -typedef ^ ParameterType IntKi DirIndexDeltaHigh - - - "Directory index delta for high-resolution AMReX wind" +typedef ^ ParameterType IntKi DirIndexLow {:} - - "AMReX directory index for each low-resolution time step; 0-based, indexed by the time step number n" - +typedef ^ ParameterType IntKi DirIndexHigh {:} - - "AMReX directory index for each high-resolution time step; 0-based, shared by all high-resolution sub-volumes" - typedef ^ ParameterType InflowWind_ParameterType IfW {:} - - "InflowWind module parameters" - # parameters for output #typedef ^ ParameterType IntKi NumOuts - - - "Number of parameters in the output list (number of outputs requested)" - @@ -225,8 +227,12 @@ typedef ^ ParameterType IntKi NOutDisWindXZ - - - "Number typedef ^ ParameterType ReKi OutDisWindY {:} - - "Y coordinates of XZ planes for output of disturbed wind data across the low-resolution domain [1 to NOutDisWindXZ]" meters typedef ^ ParameterType LOGICAL OutDisWindYvalid {:} - - "Valid XZ planes for output of disturbed wind data across the low-resolution domain [1 to NOutDisWindXZ]" - typedef ^ ParameterType CHARACTER(1024) OutFileRoot - - - "The root name derived from the primary FAST.Farm input file" - -typedef ^ ParameterType CHARACTER(1024) OutFileVTKRoot - - - "The root name for VTK outputs" - +typedef ^ ParameterType CHARACTER(1024) OutFileFFvtkRoot - - - "The root name for VTK outputs" - +typedef ^ ParameterType CHARACTER(1024) OutFileFFvtkWakeRoot - - - "The root name for VTK outputs for wake planes" - +typedef ^ ParameterType CHARACTER(1024) OutFileFFvtkWakeNullData - - - "Null data for an unpopulated wake plane" - typedef ^ ParameterType IntKi VTK_tWidth - - - "Number of characters for VTK timestamp outputs" - +typedef ^ ParameterType IntKi VTK_tWidthPlanes - 0 - "Number of charactes for the VTK plane numbers" - +typedef ^ ParameterType LOGICAL WrPlanes - .false. - "Write plane data out" - #typedef ^ ParameterType OutParmType OutParam {:} - - "Names and units (and other characteristics) of all requested output parameters" - #wake added turbulence typedef ^ ParameterType Logical WAT_Enabled - - - "Switch for turning on and off wake-added turbulence" - diff --git a/modules/awae/src/AWAE_Types.f90 b/modules/awae/src/AWAE_Types.f90 index ad08ebfd35..57a88f6afc 100644 --- a/modules/awae/src/AWAE_Types.f90 +++ b/modules/awae/src/AWAE_Types.f90 @@ -65,6 +65,7 @@ MODULE AWAE_Types INTEGER(IntKi) :: NOutDisWindXZ = 0_IntKi !< Number of XZ planes for output of disturbed wind data across the low-resolution domain to /Low/DisXZ..t.vtk [0 to 9] [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: OutDisWindY !< Y coordinates of XZ planes for output of disturbed wind data across the low-resolution domain [1 to NOutDisWindXZ] [meters] REAL(DbKi) :: WrDisDT = 0.0_R8Ki !< The time between vtk outputs [must be a multiple of the low resolution time step] [s] + LOGICAL :: OutAllPlanes = .false. !< Output all wake planes in VTK at all time steps [-] LOGICAL :: ChkWndFiles = .false. !< Check all the ambient wind files for data consistency (flag) [-] INTEGER(IntKi) :: Mod_Meander = 0_IntKi !< Spatial filter model for wake meandering {1: uniform, 2: truncated jinc, 3: windowed jinc} [DEFAULT=2] [-] REAL(ReKi) :: C_Meander = 0.0_ReKi !< Calibrated parameter for wake meandering [>=1.0] [DEFAULT=1.9] [-] @@ -154,13 +155,13 @@ MODULE AWAE_Types INTEGER(IntKi) , DIMENSION(:,:,:), ALLOCATABLE :: iPlaneTurbChunk !< First and Last plane index by source turbine and destination chunk index [-] LOGICAL , DIMENSION(:), ALLOCATABLE :: LowResChunkHasWake !< Low-res gridFirst and Last plane index by source turbine and destination chunk index [-] REAL(ReKi) :: MaxWakePointSep = 0.0_ReKi !< Maximum separation between wake points [-] - LOGICAL , DIMENSION(:,:), ALLOCATABLE :: parallelFlag !< [-] - REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: r_s !< [-] - REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: r_e !< [-] - REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: rhat_s !< [-] - REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: rhat_e !< [-] - REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: pvec_cs !< [-] - REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: pvec_ce !< [-] + LOGICAL , DIMENSION(:,:), ALLOCATABLE :: parallelFlag !< Flag indicating whether adjacent wake planes np and np+1 are (numerically) parallel for each turbine; dims: (plane-pair index np, turbine index nt) [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: r_s !< Perpendicular distance from the start wake-plane center p_plane(:,np,nt) to the intersection line of adjacent (non-parallel) wake planes np and np+1; dims: (plane-pair index np, turbine index nt) [m] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: r_e !< Perpendicular distance from the end wake-plane center p_plane(:,np+1,nt) to the intersection line of adjacent (non-parallel) wake planes np and np+1; dims: (plane-pair index np, turbine index nt) [m] + REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: rhat_s !< Unit vector lying in the start wake plane (np) pointing from the plane-plane intersection line toward the start plane center; dims: (XYZ component, plane-pair index np, turbine index nt) [-] + REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: rhat_e !< Unit vector lying in the end wake plane (np+1) pointing from the plane-plane intersection line toward the end plane center; dims: (XYZ component, plane-pair index np, turbine index nt) [-] + REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: pvec_cs !< Closest point on the plane-plane intersection line to the start wake-plane center, p_plane(:,np,nt) - r_s*rhat_s; dims: (XYZ component, plane-pair index np, turbine index nt) [m] + REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: pvec_ce !< Closest point on the plane-plane intersection line to the end wake-plane center, p_plane(:,np+1,nt) - r_e*rhat_e; dims: (XYZ component, plane-pair index np, turbine index nt) [m] REAL(SiKi) , DIMENSION(:,:,:,:), ALLOCATABLE :: outVizXYPlane !< An array holding the output data for a 2D visualization slice [-] REAL(SiKi) , DIMENSION(:,:,:,:), ALLOCATABLE :: outVizYZPlane !< An array holding the output data for a 2D visualization slice [-] REAL(SiKi) , DIMENSION(:,:,:,:), ALLOCATABLE :: outVizXZPlane !< An array holding the output data for a 2D visualization slice [-] @@ -169,7 +170,8 @@ MODULE AWAE_Types TYPE(InflowWind_OutputType) :: y_IfW_Low !< InflowWind module outputs for the low-resolution grid [-] TYPE(InflowWind_OutputType) , DIMENSION(:), ALLOCATABLE :: y_IfW_High !< InflowWind module outputs for the high-resolution grid [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: V_amb_low_disk !< Rotor averaged ambiend wind speed for each wind turbine (3 x nWT) [m/s] - INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: planeDomainExit !< Value indicates edge number (0: still in domain, +/-1: +/-X, +/-2: +/-Y, +/-3: +/-Z) the plane crossed [-] + INTEGER(IntKi) , DIMENSION(:,:,:), ALLOCATABLE :: planeDomainExit !< Per-dimension flag (0: still in domain, -1: crossed lower bound, +1: crossed upper bound) for each plane [dim,plane,turbine] [-] + INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: WakeVTK_StartN !< VTK output index (n/WrDisSkp1) when wake plane is first written. Indices [wakenum,turbnum] [-] END TYPE AWAE_MiscVarType ! ======================= ! ========= LRGChunkType ======= @@ -194,7 +196,6 @@ MODULE AWAE_Types INTEGER(IntKi) :: nPoints = 0_IntKi !< Number of spatial nodes [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: GridPoints !< XYZ components (global positions) of the spatial discretization of the grid [m] REAL(ReKi) , DIMENSION(1:3) :: Size = 0.0_ReKi !< XYZ size of the grid [m] - REAL(ReKi) , DIMENSION(1:3) :: Center = 0.0_ReKi !< XYZ coordinates of the grid center [m] TYPE(LRGChunkType) , DIMENSION(:), ALLOCATABLE :: WakeChunks !< Chunks for updating grid from wake [-] END TYPE LRGParamType ! ======================= @@ -235,9 +236,8 @@ MODULE AWAE_Types INTEGER(IntKi) :: Mod_Projection = 0_IntKi !< Switch to select how the wake plane velocity is projected in AWAE {1: keep all components, 2: project against plane normal} or DEFAULT [DEFAULT=1: if Mod_Wake is 1 or 3, or DEFAULT=2: if Mod_Wake is 2] [-] character(12) :: DirStartIndex !< Starting directory index suffix for AMReX wind [-] INTEGER(IntKi) :: DirIndexLen = 0_IntKi !< Number of characters in directory index [-] - INTEGER(IntKi) :: DirStartNum = 0_IntKi !< Starting directory index number for AMReX wind [-] - INTEGER(IntKi) :: DirIndexDeltaLow = 0_IntKi !< Directory index delta for low-resolution AMReX wind [-] - INTEGER(IntKi) :: DirIndexDeltaHigh = 0_IntKi !< Directory index delta for high-resolution AMReX wind [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: DirIndexLow !< AMReX directory index for each low-resolution time step; 0-based, indexed by the time step number n [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: DirIndexHigh !< AMReX directory index for each high-resolution time step; 0-based, shared by all high-resolution sub-volumes [-] TYPE(InflowWind_ParameterType) , DIMENSION(:), ALLOCATABLE :: IfW !< InflowWind module parameters [-] INTEGER(IntKi) :: WrDisSkp1 = 0_IntKi !< Number of time steps to skip plus one [-] LOGICAL :: WrDisWind = .false. !< Write disturbed wind data to /Low/Dis.t.vtk etc.? [-] @@ -251,8 +251,12 @@ MODULE AWAE_Types REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: OutDisWindY !< Y coordinates of XZ planes for output of disturbed wind data across the low-resolution domain [1 to NOutDisWindXZ] [meters] LOGICAL , DIMENSION(:), ALLOCATABLE :: OutDisWindYvalid !< Valid XZ planes for output of disturbed wind data across the low-resolution domain [1 to NOutDisWindXZ] [-] CHARACTER(1024) :: OutFileRoot !< The root name derived from the primary FAST.Farm input file [-] - CHARACTER(1024) :: OutFileVTKRoot !< The root name for VTK outputs [-] + CHARACTER(1024) :: OutFileFFvtkRoot !< The root name for VTK outputs [-] + CHARACTER(1024) :: OutFileFFvtkWakeRoot !< The root name for VTK outputs for wake planes [-] + CHARACTER(1024) :: OutFileFFvtkWakeNullData !< Null data for an unpopulated wake plane [-] INTEGER(IntKi) :: VTK_tWidth = 0_IntKi !< Number of characters for VTK timestamp outputs [-] + INTEGER(IntKi) :: VTK_tWidthPlanes = 0 !< Number of charactes for the VTK plane numbers [-] + LOGICAL :: WrPlanes = .false. !< Write plane data out [-] LOGICAL :: WAT_Enabled = .false. !< Switch for turning on and off wake-added turbulence [-] TYPE(FlowFieldType) , POINTER :: WAT_FlowField => NULL() !< Pointer to the InflowWinds flow field data type [-] END TYPE AWAE_ParameterType @@ -458,6 +462,7 @@ subroutine AWAE_CopyInputFileType(SrcInputFileTypeData, DstInputFileTypeData, Ct DstInputFileTypeData%OutDisWindY = SrcInputFileTypeData%OutDisWindY end if DstInputFileTypeData%WrDisDT = SrcInputFileTypeData%WrDisDT + DstInputFileTypeData%OutAllPlanes = SrcInputFileTypeData%OutAllPlanes DstInputFileTypeData%ChkWndFiles = SrcInputFileTypeData%ChkWndFiles DstInputFileTypeData%Mod_Meander = SrcInputFileTypeData%Mod_Meander DstInputFileTypeData%C_Meander = SrcInputFileTypeData%C_Meander @@ -621,6 +626,7 @@ subroutine AWAE_PackInputFileType(RF, Indata) call RegPack(RF, InData%NOutDisWindXZ) call RegPackAlloc(RF, InData%OutDisWindY) call RegPack(RF, InData%WrDisDT) + call RegPack(RF, InData%OutAllPlanes) call RegPack(RF, InData%ChkWndFiles) call RegPack(RF, InData%Mod_Meander) call RegPack(RF, InData%C_Meander) @@ -672,6 +678,7 @@ subroutine AWAE_UnPackInputFileType(RF, OutData) call RegUnpack(RF, OutData%NOutDisWindXZ); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%OutDisWindY); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WrDisDT); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%OutAllPlanes); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%ChkWndFiles); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Mod_Meander); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%C_Meander); if (RegCheckErr(RF, RoutineName)) return @@ -1441,10 +1448,10 @@ subroutine AWAE_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) DstMiscData%V_amb_low_disk = SrcMiscData%V_amb_low_disk end if if (allocated(SrcMiscData%planeDomainExit)) then - LB(1:2) = lbound(SrcMiscData%planeDomainExit) - UB(1:2) = ubound(SrcMiscData%planeDomainExit) + LB(1:3) = lbound(SrcMiscData%planeDomainExit) + UB(1:3) = ubound(SrcMiscData%planeDomainExit) if (.not. allocated(DstMiscData%planeDomainExit)) then - allocate(DstMiscData%planeDomainExit(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + allocate(DstMiscData%planeDomainExit(LB(1):UB(1),LB(2):UB(2),LB(3):UB(3)), stat=ErrStat2) if (ErrStat2 /= 0) then call SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%planeDomainExit.', ErrStat, ErrMsg, RoutineName) return @@ -1452,6 +1459,18 @@ subroutine AWAE_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) end if DstMiscData%planeDomainExit = SrcMiscData%planeDomainExit end if + if (allocated(SrcMiscData%WakeVTK_StartN)) then + LB(1:2) = lbound(SrcMiscData%WakeVTK_StartN) + UB(1:2) = ubound(SrcMiscData%WakeVTK_StartN) + if (.not. allocated(DstMiscData%WakeVTK_StartN)) then + allocate(DstMiscData%WakeVTK_StartN(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%WakeVTK_StartN.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstMiscData%WakeVTK_StartN = SrcMiscData%WakeVTK_StartN + end if end subroutine subroutine AWAE_DestroyMisc(MiscData, ErrStat, ErrMsg) @@ -1564,6 +1583,9 @@ subroutine AWAE_DestroyMisc(MiscData, ErrStat, ErrMsg) if (allocated(MiscData%planeDomainExit)) then deallocate(MiscData%planeDomainExit) end if + if (allocated(MiscData%WakeVTK_StartN)) then + deallocate(MiscData%WakeVTK_StartN) + end if end subroutine subroutine AWAE_PackMisc(RF, Indata) @@ -1626,6 +1648,7 @@ subroutine AWAE_PackMisc(RF, Indata) end if call RegPackAlloc(RF, InData%V_amb_low_disk) call RegPackAlloc(RF, InData%planeDomainExit) + call RegPackAlloc(RF, InData%WakeVTK_StartN) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -1703,6 +1726,7 @@ subroutine AWAE_UnPackMisc(RF, OutData) end if call RegUnpackAlloc(RF, OutData%V_amb_low_disk); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%planeDomainExit); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%WakeVTK_StartN); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine AWAE_CopyLRGChunkType(SrcLRGChunkTypeData, DstLRGChunkTypeData, CtrlCode, ErrStat, ErrMsg) @@ -1819,7 +1843,6 @@ subroutine AWAE_CopyLRGParamType(SrcLRGParamTypeData, DstLRGParamTypeData, CtrlC DstLRGParamTypeData%GridPoints = SrcLRGParamTypeData%GridPoints end if DstLRGParamTypeData%Size = SrcLRGParamTypeData%Size - DstLRGParamTypeData%Center = SrcLRGParamTypeData%Center if (allocated(SrcLRGParamTypeData%WakeChunks)) then LB(1:1) = lbound(SrcLRGParamTypeData%WakeChunks) UB(1:1) = ubound(SrcLRGParamTypeData%WakeChunks) @@ -1876,7 +1899,6 @@ subroutine AWAE_PackLRGParamType(RF, Indata) call RegPack(RF, InData%nPoints) call RegPackAlloc(RF, InData%GridPoints) call RegPack(RF, InData%Size) - call RegPack(RF, InData%Center) call RegPack(RF, allocated(InData%WakeChunks)) if (allocated(InData%WakeChunks)) then call RegPackBounds(RF, 1, lbound(InData%WakeChunks), ubound(InData%WakeChunks)) @@ -1904,7 +1926,6 @@ subroutine AWAE_UnPackLRGParamType(RF, OutData) call RegUnpack(RF, OutData%nPoints); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%GridPoints); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Size); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%Center); if (RegCheckErr(RF, RoutineName)) return if (allocated(OutData%WakeChunks)) deallocate(OutData%WakeChunks) call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return if (IsAllocAssoc) then @@ -2075,9 +2096,30 @@ subroutine AWAE_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%Mod_Projection = SrcParamData%Mod_Projection DstParamData%DirStartIndex = SrcParamData%DirStartIndex DstParamData%DirIndexLen = SrcParamData%DirIndexLen - DstParamData%DirStartNum = SrcParamData%DirStartNum - DstParamData%DirIndexDeltaLow = SrcParamData%DirIndexDeltaLow - DstParamData%DirIndexDeltaHigh = SrcParamData%DirIndexDeltaHigh + if (allocated(SrcParamData%DirIndexLow)) then + LB(1:1) = lbound(SrcParamData%DirIndexLow) + UB(1:1) = ubound(SrcParamData%DirIndexLow) + if (.not. allocated(DstParamData%DirIndexLow)) then + allocate(DstParamData%DirIndexLow(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%DirIndexLow.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%DirIndexLow = SrcParamData%DirIndexLow + end if + if (allocated(SrcParamData%DirIndexHigh)) then + LB(1:1) = lbound(SrcParamData%DirIndexHigh) + UB(1:1) = ubound(SrcParamData%DirIndexHigh) + if (.not. allocated(DstParamData%DirIndexHigh)) then + allocate(DstParamData%DirIndexHigh(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%DirIndexHigh.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%DirIndexHigh = SrcParamData%DirIndexHigh + end if if (allocated(SrcParamData%IfW)) then LB(1:1) = lbound(SrcParamData%IfW) UB(1:1) = ubound(SrcParamData%IfW) @@ -2172,8 +2214,12 @@ subroutine AWAE_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%OutDisWindYvalid = SrcParamData%OutDisWindYvalid end if DstParamData%OutFileRoot = SrcParamData%OutFileRoot - DstParamData%OutFileVTKRoot = SrcParamData%OutFileVTKRoot + DstParamData%OutFileFFvtkRoot = SrcParamData%OutFileFFvtkRoot + DstParamData%OutFileFFvtkWakeRoot = SrcParamData%OutFileFFvtkWakeRoot + DstParamData%OutFileFFvtkWakeNullData = SrcParamData%OutFileFFvtkWakeNullData DstParamData%VTK_tWidth = SrcParamData%VTK_tWidth + DstParamData%VTK_tWidthPlanes = SrcParamData%VTK_tWidthPlanes + DstParamData%WrPlanes = SrcParamData%WrPlanes DstParamData%WAT_Enabled = SrcParamData%WAT_Enabled DstParamData%WAT_FlowField => SrcParamData%WAT_FlowField end subroutine @@ -2206,6 +2252,12 @@ subroutine AWAE_DestroyParam(ParamData, ErrStat, ErrMsg) if (allocated(ParamData%z)) then deallocate(ParamData%z) end if + if (allocated(ParamData%DirIndexLow)) then + deallocate(ParamData%DirIndexLow) + end if + if (allocated(ParamData%DirIndexHigh)) then + deallocate(ParamData%DirIndexHigh) + end if if (allocated(ParamData%IfW)) then LB(1:1) = lbound(ParamData%IfW) UB(1:1) = ubound(ParamData%IfW) @@ -2274,9 +2326,8 @@ subroutine AWAE_PackParam(RF, Indata) call RegPack(RF, InData%Mod_Projection) call RegPack(RF, InData%DirStartIndex) call RegPack(RF, InData%DirIndexLen) - call RegPack(RF, InData%DirStartNum) - call RegPack(RF, InData%DirIndexDeltaLow) - call RegPack(RF, InData%DirIndexDeltaHigh) + call RegPackAlloc(RF, InData%DirIndexLow) + call RegPackAlloc(RF, InData%DirIndexHigh) call RegPack(RF, allocated(InData%IfW)) if (allocated(InData%IfW)) then call RegPackBounds(RF, 1, lbound(InData%IfW), ubound(InData%IfW)) @@ -2298,8 +2349,12 @@ subroutine AWAE_PackParam(RF, Indata) call RegPackAlloc(RF, InData%OutDisWindY) call RegPackAlloc(RF, InData%OutDisWindYvalid) call RegPack(RF, InData%OutFileRoot) - call RegPack(RF, InData%OutFileVTKRoot) + call RegPack(RF, InData%OutFileFFvtkRoot) + call RegPack(RF, InData%OutFileFFvtkWakeRoot) + call RegPack(RF, InData%OutFileFFvtkWakeNullData) call RegPack(RF, InData%VTK_tWidth) + call RegPack(RF, InData%VTK_tWidthPlanes) + call RegPack(RF, InData%WrPlanes) call RegPack(RF, InData%WAT_Enabled) call RegPack(RF, associated(InData%WAT_FlowField)) if (associated(InData%WAT_FlowField)) then @@ -2356,9 +2411,8 @@ subroutine AWAE_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%Mod_Projection); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%DirStartIndex); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%DirIndexLen); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%DirStartNum); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%DirIndexDeltaLow); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%DirIndexDeltaHigh); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%DirIndexLow); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%DirIndexHigh); if (RegCheckErr(RF, RoutineName)) return if (allocated(OutData%IfW)) deallocate(OutData%IfW) call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return if (IsAllocAssoc) then @@ -2384,8 +2438,12 @@ subroutine AWAE_UnPackParam(RF, OutData) call RegUnpackAlloc(RF, OutData%OutDisWindY); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%OutDisWindYvalid); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%OutFileRoot); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%OutFileVTKRoot); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%OutFileFFvtkRoot); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%OutFileFFvtkWakeRoot); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%OutFileFFvtkWakeNullData); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%VTK_tWidth); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%VTK_tWidthPlanes); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%WrPlanes); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WAT_Enabled); if (RegCheckErr(RF, RoutineName)) return if (associated(OutData%WAT_FlowField)) deallocate(OutData%WAT_FlowField) call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return diff --git a/modules/awae/src/AWAE_vtk.f90 b/modules/awae/src/AWAE_vtk.f90 new file mode 100644 index 0000000000..4ff1ca531c --- /dev/null +++ b/modules/awae/src/AWAE_vtk.f90 @@ -0,0 +1,509 @@ +!********************************************************************************************************************************** +! LICENSING +! Copyright (C) 2015-2016 National Renewable Energy Laboratory +! +! This file is part of Ambient Wind and Array Effects model for FAST.Farm. +! +! Licensed under the Apache License, Version 2.0 (the "License"); +! you may not use this file except in compliance with the License. +! You may obtain a copy of the License at +! +! http://www.apache.org/licenses/LICENSE-2.0 +! +! Unless required by applicable law or agreed to in writing, software +! distributed under the License is distributed on an "AS IS" BASIS, +! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +! See the License for the specific language governing permissions and +! limitations under the License. +! +!********************************************************************************************************************************** +!> VTK output routines for the AWAE module (wake-plane wireframes, structured +!! grid data files, and ParaView .vtk.series index files). +module AWAE_vtk + + use NWTC_Library + use AWAE_Types + use VTK + +#ifdef _OPENMP + use OMP_LIB +#endif + + implicit none + + private + + public :: PlaneAxes + public :: PlaneCorners + public :: Write_Planes_WireFrame + public :: Write_Planes_Data + public :: Write_WakePlane_Data_File + public :: Write_WakePlane_Series + public :: Write_WireFrame_Series + public :: Write_DisWind_Series + public :: Write_NullPlane + +contains + +!> Construct the orthonormal in-plane basis (yhat, zhat) in the global +!! inertial frame for a wake plane with unit normal `xhat`. yhat is the +!! horizontal in-plane direction (no Z component); zhat = xhat x yhat. +!! If xhat is purely vertical, yhat falls back to the global Y axis. +subroutine PlaneAxes(xhat, yhat, zhat) + real(ReKi), intent(in ) :: xhat(3) + real(ReKi), intent( out) :: yhat(3) + real(ReKi), intent( out) :: zhat(3) + real(ReKi) :: ynorm + + yhat = (/ -xhat(2), xhat(1), 0.0_ReKi /) + ynorm = TwoNorm(yhat) + if (ynorm > 0.0_ReKi) then + yhat = yhat / ynorm + else + yhat = (/ 0.0_ReKi, 1.0_ReKi, 0.0_ReKi /) + end if + + zhat(1) = xhat(2)*yhat(3) - xhat(3)*yhat(2) + zhat(2) = xhat(3)*yhat(1) - xhat(1)*yhat(3) + zhat(3) = xhat(1)*yhat(2) - xhat(2)*yhat(1) +end subroutine PlaneAxes + + +!> Compute the four corners of a single wake plane in the global inertial +!! frame from the plane normal `xhat`, the plane center `pc`, and the +!! in-plane half-extents `hy` (horizontal) and `hz` (vertical-ish). +!! Corners are returned counter-clockwise about +xhat. +subroutine PlaneCorners(xhat, pc, hy, hz, corners) + real(ReKi), intent(in ) :: xhat(3) !< Plane normal (unit vector) + real(ReKi), intent(in ) :: pc(3) !< Plane center, global frame + real(ReKi), intent(in ) :: hy !< In-plane horizontal half-extent + real(ReKi), intent(in ) :: hz !< In-plane vertical-ish half-extent + real(ReKi), intent( out) :: corners(3,4) !< Four corner positions, global frame + + real(ReKi) :: yhat(3), zhat(3) + + call PlaneAxes(xhat, yhat, zhat) + + ! Four corners, ordered counter-clockwise about +xhat + corners(:,1) = pc - hy*yhat - hz*zhat + corners(:,2) = pc + hy*yhat - hz*zhat + corners(:,3) = pc + hy*yhat + hz*zhat + corners(:,4) = pc - hy*yhat + hz*zhat +end subroutine PlaneCorners + + +!> Write the four edges of every active wake plane for each turbine as a +!! wireframe to a per-turbine VTK polydata file. One file is produced per +!! turbine per output step, named .T.WakePlanesWireFrame..vtk. +subroutine Write_Planes_WireFrame(p, u, t, Tstr) + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + type(AWAE_InputType), intent(in) :: u !< AWAE inputs + real(DbKi), intent(in) :: t !< Current simulation time in seconds + character(*), intent(in) :: Tstr !< Zero-padded time-step string for filenames + + integer(IntKi) :: nt_wp, np_wp, nActive, ipt + real(ReKi) :: pc(3), corners(3,4) + real(ReKi) :: hy, hz + real(ReKi), allocatable :: WPPoints(:,:) + integer(IntKi), allocatable :: WPLines(:,:) + type(VTK_Misc) :: mvtk + character(1024) :: WPFileName + integer(IntKi) :: ntWidth + character(16) :: TurbNum, FmtStrT + + ! Plane half-extents in the local Y and Z (in-plane) directions + hy = p%y(p%NumRadii-1) + hz = p%z(p%NumRadii-1) + + ! Zero-padded width sufficient to represent all turbine indices + ntWidth = max(1, int(floor(log10(real(max(p%NumTurbines, 1), ReKi))) + 1, IntKi)) + write(FmtStrT, '(A,I0,A,I0,A)') '(I', ntWidth, '.', ntWidth, ')' + + do nt_wp = 1, p%NumTurbines + nActive = NINT(u%NumPlanes(nt_wp)) + if (nActive <= 0) cycle + + allocate(WPPoints(3, 4*nActive)) + allocate(WPLines(2, 4*nActive)) + + do np_wp = 0, nActive - 1 + pc = u%p_plane(:, np_wp, nt_wp) + call PlaneCorners(u%xhat_plane(:, np_wp, nt_wp), pc, hy, hz, corners) + + ipt = 4*np_wp + WPPoints(:, ipt+1) = corners(:,1) + WPPoints(:, ipt+2) = corners(:,2) + WPPoints(:, ipt+3) = corners(:,3) + WPPoints(:, ipt+4) = corners(:,4) + + ! Four edges of the closed quad outline (0-based point indices for VTK) + WPLines(:, ipt+1) = (/ ipt, ipt+1 /) + WPLines(:, ipt+2) = (/ ipt+1, ipt+2 /) + WPLines(:, ipt+3) = (/ ipt+2, ipt+3 /) + WPLines(:, ipt+4) = (/ ipt+3, ipt /) + end do + + write(TurbNum, FmtStrT) nt_wp + WPFileName = trim(p%OutFileFFvtkWakeRoot)//".T"//trim(TurbNum)// & + ".WakePlanesWireFrame."//trim(Tstr)//".vtk" + + call vtk_misc_init(mvtk) + if (vtk_new_ascii_file(WPFileName, & + "Wake plane wireframes for turbine "//trim(TurbNum)// & + " at time = "//trim(num2lstr(t))//" seconds.", mvtk)) then + call vtk_dataset_polydata(WPPoints, mvtk, .false.) + call vtk_lines(WPLines, mvtk) + call vtk_close_file(mvtk) + end if + + deallocate(WPPoints, WPLines) + + end do +end subroutine Write_Planes_WireFrame + + +!> Write one VTK STRUCTURED_GRID file per active wake plane containing the +!! wake velocity sampled on the plane's structured Y-Z grid. Point +!! coordinates and velocity vectors are written in the global inertial frame. +subroutine Write_Planes_Data(p, u, m, n, t, Tstr) + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + type(AWAE_InputType), intent(in) :: u !< AWAE inputs + type(AWAE_MiscVarType), intent(inout) :: m !< AWAE misc variables + integer(IntKi), intent(in) :: n !< Current low-resolution time step index + real(DbKi), intent(in) :: t !< Current simulation time in seconds + character(*), intent(in) :: Tstr !< Zero-padded time-step string for filenames + + integer(IntKi) :: nt_wp, np_wp, nActive, ntWidth + integer(IntKi) :: jy, kz, idx, nY, nZ + real(ReKi) :: pc(3), xhat(3), yhat(3), zhat(3) + real(ReKi), allocatable :: Pts(:,:), Vel(:,:) + character(1024) :: WPFileName + character(16) :: FmtStrWk, FmtStrT + character(p%VTK_tWidthPlanes) :: PlaneNum + character(6) :: TurbNumStr + + ! Plane grid dimensions (p%y and p%z run from -NumRadii+1 to NumRadii-1) + nY = 2*p%NumRadii - 1 + nZ = 2*p%NumRadii - 1 + + ! Field width for plane index, based on MaxPlanes (1-based count) + write(FmtStrWk, '(A,I0,A,I0,A)') '(I', p%VTK_tWidthPlanes, '.', p%VTK_tWidthPlanes, ')' + ! Zero-padded width sufficient to represent all turbine indices + ntWidth = max(1, int(floor(log10(real(max(p%NumTurbines, 1), ReKi))) + 1, IntKi)) + write(FmtStrT, '(A,I0,A,I0,A)') '(A1,I', ntWidth, '.', ntWidth, ')' + + allocate(Pts(3, nY*nZ)) + allocate(Vel(3, nY*nZ)) + + do nt_wp = 1, p%NumTurbines + ! Turbine number + write(TurbNumStr, FmtStrT) "T", nt_wp + + ! Number of planes active + nActive = NINT(u%NumPlanes(nt_wp)) + + do np_wp = 0, nActive - 1 + pc = u%p_plane(:, np_wp, nt_wp) + xhat = u%xhat_plane(:, np_wp, nt_wp) + call PlaneAxes(xhat, yhat, zhat) + + ! Fill points and velocity in VTK natural order (jy fastest, then kz) + do kz = 1, nZ + do jy = 1, nY + idx = jy + (kz-1)*nY + Pts(:, idx) = pc + p%y(jy-p%NumRadii) * yhat & + + p%z(kz-p%NumRadii) * zhat + Vel(:, idx) = u%Vx_wake(jy-p%NumRadii, kz-p%NumRadii, np_wp, nt_wp) * xhat & + + u%Vy_wake(jy-p%NumRadii, kz-p%NumRadii, np_wp, nt_wp) * yhat & + + u%Vz_wake(jy-p%NumRadii, kz-p%NumRadii, np_wp, nt_wp) * zhat + end do + end do + + ! Per-plane index string with consistent zero-padded width + write(PlaneNum, FmtStrWk) np_wp + + WPFileName = trim(p%OutFileFFvtkWakeRoot)//"."//trim(TurbNumStr)//".WakePlane_"//trim(PlaneNum)// & + "."//trim(Tstr)//".vtk" + + call Write_WakePlane_Data_File(WPFileName, & + "Wake plane "//trim(PlaneNum)//" at time = "// & + trim(num2lstr(t))//" seconds.", nY, nZ, Pts, Vel) + + ! track the VTK output index this plane was first written at (initialized to huge, so this will catch only the first). + ! Stored as the output index (n/WrDisSkp1), matching Tstr and the series-file loop in Write_WakePlane_Series. + if (m%WakeVTK_StartN(np_wp,nt_wp) > n/p%WrDisSkp1) m%WakeVTK_StartN(np_wp,nt_wp) = n/p%WrDisSkp1 + end do + end do + + deallocate(Pts, Vel) +end subroutine Write_Planes_Data + + +!> Helper: write a single 2D VTK STRUCTURED_GRID file for one wake plane +!! containing point coordinates and a "WakeVelocity" point-data vector. +subroutine Write_WakePlane_Data_File(WPFileName, label, n1, n2, Pts, Vel) + character(*), intent(in) :: WPFileName + character(*), intent(in) :: label + integer(IntKi), intent(in) :: n1, n2 + real(ReKi), intent(in) :: Pts(:,:) !< 3 x (n1*n2) + real(ReKi), intent(in) :: Vel(:,:) !< 3 x (n1*n2) + type(VTK_Misc) :: mvtk + + call vtk_misc_init(mvtk) + if (vtk_new_ascii_file(WPFileName, label, mvtk)) then + call vtk_dataset_structured_grid(Pts, n1, n2, 1, mvtk) + call vtk_point_data_init(mvtk) + call vtk_point_data_vector(Vel, "WakeVelocityDeficit", mvtk) + call vtk_close_file(mvtk) + end if +end subroutine Write_WakePlane_Data_File + + +!> Write the final ParaView .vtk.series index files for all wake planes +!! that were output during the simulation. Called once from AWAE_End. +!! Each series file lists every VTK output step, referencing either the +!! actual per-plane VTK file or a null-data placeholder for steps before +!! the plane first appeared. +subroutine Write_WakePlane_Series(p, m) + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + type(AWAE_MiscVarType), intent(in) :: m !< AWAE misc variables + + integer(IntKi) :: nt_wp, np_wp, ntWidth, n_final, n_out + integer(IntKi) :: UnSer, out_idx, SerErrStat + character(ErrMsgLen) :: SerErrMsg + character(16) :: FmtStrWk, FmtStrT + character(p%VTK_tWidthPlanes) :: PlaneNum + character(6) :: TurbNumStr + character(1024) :: SeriesFile, baseName, EntryName, VTKprefix + character(p%VTK_tWidth) :: TstrOut + character(32) :: TimeStr + real(DbKi) :: t_out + logical :: firstEntry + + if (.not. p%WrPlanes) return + + ! Final low-res time step and total number of VTK output steps + n_final = p%NumDT - 1 + n_out = n_final / p%WrDisSkp1 + + ! Compute basename once (shared across all series files) + call GetPath(p%OutFileFFvtkWakeRoot, EntryName, baseName) + + ! Format strings for zero-padded plane and turbine numbers + write(FmtStrWk, '(A,I0,A,I0,A)') '(I', p%VTK_tWidthPlanes, '.', p%VTK_tWidthPlanes, ')' + ntWidth = max(1, int(floor(log10(real(max(p%NumTurbines, 1), ReKi))) + 1, IntKi)) + write(FmtStrT, '(A,I0,A,I0,A)') '(A1,I', ntWidth, '.', ntWidth, ')' + + do nt_wp = 1, p%NumTurbines + write(TurbNumStr, FmtStrT) "T", nt_wp + do np_wp = 0, p%MaxPlanes - 1 + ! Skip planes that were never written during the simulation (WakeVTK_StartN is a VTK output index) + if (m%WakeVTK_StartN(np_wp, nt_wp) > n_out) cycle + + write(PlaneNum, FmtStrWk) np_wp + VTKprefix = trim(TurbNumStr)//".WakePlane_"//trim(PlaneNum) + SeriesFile = trim(p%OutFileFFvtkWakeRoot)//"."//trim(VTKprefix)//".vtk.series" + + !$OMP critical(fileopen_critical) + call GetNewUnit(UnSer, SerErrStat, SerErrMsg) + call OpenFOutFile(UnSer, SeriesFile, SerErrStat, SerErrMsg) + !$OMP end critical(fileopen_critical) + if (SerErrStat >= AbortErrLev) cycle + + write(UnSer, '(A)') '{' + write(UnSer, '(A)') ' "file-series-version" : "1.0",' + write(UnSer, '(A)') ' "files" : [' + + firstEntry = .true. + do out_idx = 0, n_out + t_out = real(out_idx, DbKi) * real(p%WrDisSkp1, DbKi) * p%dt_low + write(TstrOut, '(i'//trim(Num2LStr(p%VTK_tWidth))//'.'// & + trim(Num2LStr(p%VTK_tWidth))//')') out_idx + + if (m%WakeVTK_StartN(np_wp, nt_wp) > out_idx) then + EntryName = trim(p%OutFileFFvtkWakeNullData) + else + EntryName = trim(baseName)//"."//trim(VTKprefix)//"."//trim(TstrOut)//".vtk" + endif + + write(TimeStr, '(F14.5)') t_out + if (firstEntry) then + write(UnSer, '(A,A,A,A,A)') ' { "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + firstEntry = .false. + else + write(UnSer, '(A,A,A,A,A)') ' ,{ "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + end if + end do + + write(UnSer, '(A)') ' ]' + write(UnSer, '(A)') '}' + + !$OMP critical(fileopen_critical) + close(UnSer) + !$OMP end critical(fileopen_critical) + end do + end do +end subroutine Write_WakePlane_Series + +!> Write a zeroed-out "null" wake-plane VTK file used as a placeholder +!! in the ParaView series for time steps before a plane first appears. +subroutine Write_NullPlane(OutFileVTKwakeDir, p) + character(*), intent(in) :: OutFileVTKwakeDir !< Directory for VTK wake output files + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + + real(ReKi), allocatable :: PtsVel(:,:) + integer(IntKi) :: nY, nZ + + nY = 2*p%NumRadii - 1 + nZ = 2*p%NumRadii - 1 + allocate(PtsVel(3, nY*nZ)) + PtsVel = 0.0_ReKi + call Write_WakePlane_Data_File(trim(OutFileVTKwakeDir)//PathSep//p%OutFileFFvtkWakeNullData, & + "Wake plane NULL at time = NULL seconds.", nY, nZ, PtsVel, PtsVel) + deallocate(PtsVel) +end subroutine Write_NullPlane + + +!> Write one ParaView .vtk.series index file per turbine for the wireframe +!! VTK files produced by Write_Planes_WireFrame. Called once from AWAE_End. +subroutine Write_WireFrame_Series(p) + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + + integer(IntKi) :: nt_wp, ntWidth, n_out + integer(IntKi) :: UnSer, out_idx, SerErrStat + character(ErrMsgLen) :: SerErrMsg + character(16) :: FmtStrT, TurbNum + character(1024) :: SeriesFile, baseName, EntryName + character(p%VTK_tWidth) :: TstrOut + character(32) :: TimeStr + real(DbKi) :: t_out + logical :: firstEntry + + if (.not. p%WrPlanes) return + + ! Total number of VTK output steps + n_out = (p%NumDT - 1) / p%WrDisSkp1 + + ! Compute basename once (filename portion of OutFileFFvtkWakeRoot) + call GetPath(p%OutFileFFvtkWakeRoot, EntryName, baseName) + + ! Format string for zero-padded turbine number + ntWidth = max(1, int(floor(log10(real(max(p%NumTurbines, 1), ReKi))) + 1, IntKi)) + write(FmtStrT, '(A,I0,A,I0,A)') '(I', ntWidth, '.', ntWidth, ')' + + do nt_wp = 1, p%NumTurbines + write(TurbNum, FmtStrT) nt_wp + + SeriesFile = trim(p%OutFileFFvtkWakeRoot)//".T"//trim(TurbNum)// & + ".WakePlanesWireFrame.vtk.series" + + !$OMP critical(fileopen_critical) + call GetNewUnit(UnSer, SerErrStat, SerErrMsg) + call OpenFOutFile(UnSer, SeriesFile, SerErrStat, SerErrMsg) + !$OMP end critical(fileopen_critical) + if (SerErrStat >= AbortErrLev) cycle + + write(UnSer, '(A)') '{' + write(UnSer, '(A)') ' "file-series-version" : "1.0",' + write(UnSer, '(A)') ' "files" : [' + + firstEntry = .true. + do out_idx = 0, n_out + t_out = real(out_idx, DbKi) * real(p%WrDisSkp1, DbKi) * p%dt_low + write(TstrOut, '(i'//trim(Num2LStr(p%VTK_tWidth))//'.'// & + trim(Num2LStr(p%VTK_tWidth))//')') out_idx + + EntryName = trim(baseName)//".T"//trim(TurbNum)// & + ".WakePlanesWireFrame."//trim(TstrOut)//".vtk" + + write(TimeStr, '(F14.5)') t_out + if (firstEntry) then + write(UnSer, '(A,A,A,A,A)') ' { "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + firstEntry = .false. + else + write(UnSer, '(A,A,A,A,A)') ' ,{ "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + end if + end do + + write(UnSer, '(A)') ' ]' + write(UnSer, '(A)') '}' + + !$OMP critical(fileopen_critical) + close(UnSer) + !$OMP end critical(fileopen_critical) + end do +end subroutine Write_WireFrame_Series + + +!> Write a ParaView .vtk.series index file for one disturbed-wind output +!! slice (XY, YZ, or XZ). The slice label and 1-based index select which +!! series of per-timestep VTK files to reference. +subroutine Write_DisWind_Series(p, sliceLabel, iSlice) + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + character(*), intent(in) :: sliceLabel !< "DisXY", "DisYZ", or "DisXZ" + integer(IntKi), intent(in) :: iSlice !< 1-based slice index + + integer(IntKi) :: n_out + integer(IntKi) :: UnSer, out_idx, SerErrStat + character(ErrMsgLen) :: SerErrMsg + character(1024) :: SeriesFile, baseName, EntryName + character(p%VTK_tWidth) :: TstrOut + character(32) :: TimeStr + character(3) :: PlaneNumStr + real(DbKi) :: t_out + logical :: firstEntry + + ! Total number of VTK output steps + n_out = (p%NumDT - 1) / p%WrDisSkp1 + + ! 3-digit zero-padded slice index (matches CalcOutput naming) + write(PlaneNumStr, '(i3.3)') iSlice + + ! Compute basename (filename portion of OutFileFFvtkRoot) + call GetPath(p%OutFileFFvtkRoot, EntryName, baseName) + + ! Series filename + SeriesFile = trim(p%OutFileFFvtkRoot)//".Low."//trim(sliceLabel)//PlaneNumStr//".vtk.series" + + !$OMP critical(fileopen_critical) + call GetNewUnit(UnSer, SerErrStat, SerErrMsg) + call OpenFOutFile(UnSer, SeriesFile, SerErrStat, SerErrMsg) + !$OMP end critical(fileopen_critical) + if (SerErrStat >= AbortErrLev) return + + write(UnSer, '(A)') '{' + write(UnSer, '(A)') ' "file-series-version" : "1.0",' + write(UnSer, '(A)') ' "files" : [' + + firstEntry = .true. + do out_idx = 0, n_out + t_out = real(out_idx, DbKi) * real(p%WrDisSkp1, DbKi) * p%dt_low + write(TstrOut, '(i'//trim(Num2LStr(p%VTK_tWidth))//'.'// & + trim(Num2LStr(p%VTK_tWidth))//')') out_idx + + EntryName = trim(baseName)//".Low."//trim(sliceLabel)//PlaneNumStr//"."//trim(TstrOut)//".vtk" + + write(TimeStr, '(F14.5)') t_out + if (firstEntry) then + write(UnSer, '(A,A,A,A,A)') ' { "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + firstEntry = .false. + else + write(UnSer, '(A,A,A,A,A)') ' ,{ "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + end if + end do + + write(UnSer, '(A)') ' ]' + write(UnSer, '(A)') '}' + + !$OMP critical(fileopen_critical) + close(UnSer) + !$OMP end critical(fileopen_critical) +end subroutine Write_DisWind_Series + +end module AWAE_vtk diff --git a/modules/awae/src/amrex_utils.F90 b/modules/awae/src/amrex_utils.F90 index 26d271bfe1..c66b9e816a 100644 --- a/modules/awae/src/amrex_utils.F90 +++ b/modules/awae/src/amrex_utils.F90 @@ -24,18 +24,30 @@ subroutine amrex_read_header_c(dir_path, t, dims, dx, origin, err_stat, err_msg, integer(kind=c_int), intent(in) :: err_msg_len end subroutine - subroutine amrex_read_data_c(dir_path, data, err_stat, err_msg, err_msg_len) bind(c) + subroutine amrex_read_data_c(dir_path, data, dims_expected, err_stat, err_msg, err_msg_len) bind(c) import implicit none character(kind=c_char), intent(in) :: dir_path(*) real(kind=c_float), intent(out) :: data(*) + integer(kind=c_int), intent(in) :: dims_expected(3) integer(kind=c_int), intent(out) :: err_stat character(kind=c_char), intent(out) :: err_msg(*) integer(kind=c_int), intent(in) :: err_msg_len end subroutine - subroutine amrex_find_subvols_c(dir_path, subvol, dt, num_step, start_index, first_index, & - index_delta, err_stat, err_msg, err_msg_len) bind(c) + subroutine amrex_header_text_c(dir_path, t, dims, dx, origin, ok) bind(c) + import + implicit none + character(kind=c_char), intent(in) :: dir_path(*) + real(kind=c_double), intent(out) :: t + integer(kind=c_int), intent(out) :: dims(3) + real(kind=c_double), intent(out) :: dx(3) + real(kind=c_double), intent(out) :: origin(3) + integer(kind=c_int), intent(out) :: ok + end subroutine + + subroutine amrex_find_subvols_c(dir_path, subvol, dt, num_step, start_index, dir_indices, & + err_stat, err_msg, err_msg_len) bind(c) import implicit none character(kind=c_char), intent(in) :: dir_path(*) @@ -43,8 +55,7 @@ subroutine amrex_find_subvols_c(dir_path, subvol, dt, num_step, start_index, fir character(kind=c_char), intent(in) :: start_index(*) real(kind=c_double), intent(in) :: dt integer(kind=c_int), intent(in) :: num_step - integer(kind=c_int), intent(out) :: first_index - integer(kind=c_int), intent(out) :: index_delta + integer(kind=c_int), intent(out) :: dir_indices(*) integer(kind=c_int), intent(out) :: err_stat character(kind=c_char), intent(out) :: err_msg(*) integer(kind=c_int), intent(in) :: err_msg_len @@ -52,7 +63,7 @@ subroutine amrex_find_subvols_c(dir_path, subvol, dt, num_step, start_index, fir end interface public :: amrex_init, amrex_finalize -public :: amrex_read_header, amrex_read_data, amrex_find_subvols +public :: amrex_read_header, amrex_read_data, amrex_find_subvols, amrex_parse_header_text contains @@ -84,6 +95,9 @@ subroutine amrex_read_header(DirPath, time, nXYZ, dXYZ, oXYZ, ErrStat, ErrMsg) real(c_double) :: t, origin(3), gridSpacing(3) integer(IntKi) :: i + ErrStat = ErrID_None + ErrMsg = "" + #ifdef ENABLE_AMREX_LIB ! Convert directory path to C type @@ -116,19 +130,33 @@ subroutine amrex_read_data(DirPath, gridData, ErrStat, ErrMsg) integer(IntKi), intent(out) :: ErrStat character(*), intent(out) :: ErrMsg + character(*), parameter :: RoutineName = 'amrex_read_data' character(c_char), allocatable :: dir_path(:) integer(c_int) :: err_stat_c + integer(c_int) :: dims_expected(3) character(c_char) :: err_msg_c(ErrMsgLen) integer(IntKi) :: i + ErrStat = ErrID_None + ErrMsg = "" + #ifdef ENABLE_AMREX_LIB + ! The C routine writes gridData by grid index, so it needs the extent of the destination to + ! stay inside it. The first dimension is the three velocity components. + if (size(gridData, 1) /= 3) then + call SetErrStat(ErrID_Fatal, "gridData must have 3 velocity components, got "// & + trim(Num2LStr(size(gridData, 1))), ErrStat, ErrMsg, RoutineName) + return + end if + dims_expected = int([size(gridData, 2), size(gridData, 3), size(gridData, 4)], c_int) + ! Convert directory path to C type allocate(dir_path(len_trim(DirPath) + 1)) dir_path = transfer(trim(DirPath) // c_null_char, dir_path) - ! Call C++ function to read header - call amrex_read_data_c(dir_path, gridData, err_stat_c, err_msg_c, ErrMsgLen) + ! Call C++ function to read the grid data + call amrex_read_data_c(dir_path, gridData, dims_expected, err_stat_c, err_msg_c, ErrMsgLen) ! Transfer outputs back to fortran types ErrStat = int(err_stat_c, IntKi) @@ -141,33 +169,46 @@ subroutine amrex_read_data(DirPath, gridData, ErrStat, ErrMsg) #endif end subroutine -! Search for AMReX directories based on given directory prefix, sub-volume -! number, time step, total number of steps, and the starting index string -! (e.g. `00000`). This function returns first index as a number, and the delta -! between successive directory indices. It also checks that a sufficient number -! of directories are available, that the data matches the requested time step, -! and the grid properties are consistent (size, origin, spacing). +! Search for AMReX directories based on given directory prefix, sub-volume +! number, time step, total number of steps, and the starting index string +! (e.g. `00000`). Returns DirIndices(0:NumStep-1), the directory index suffix +! to use for each time step. Directories are matched to time steps by the +! simulation time recorded in their Header; no constant stride between +! successive directory indices is assumed, so LES output written with a varying +! solver time step is supported. Also checks that every requested step is +! represented -- and, among the directories scanned, represented only once -- +! and that the grid properties are consistent (size, origin, spacing). subroutine amrex_find_subvols(DirPath, SubVol, DT, NumStep, StartIndex, & - FirstIndex, IndexDelta, ErrStat, ErrMsg) - character(*), intent(in) :: DirPath - integer(IntKi), intent(in) :: SubVol - real(DbKi), intent(in) :: DT ! Time step - integer(IntKi), intent(in) :: NumStep ! Number of steps - character(*), intent(in) :: StartIndex - integer(IntKi), intent(out) :: FirstIndex, IndexDelta - integer(IntKi), intent(out) :: ErrStat - character(*), intent(out) :: ErrMsg - + DirIndices, ErrStat, ErrMsg) + character(*), intent(in) :: DirPath + integer(IntKi), intent(in) :: SubVol + real(DbKi), intent(in) :: DT ! Time step + integer(IntKi), intent(in) :: NumStep ! Number of steps + character(*), intent(in) :: StartIndex + integer(IntKi), allocatable, intent(out) :: DirIndices(:) ! (0:NumStep-1) + integer(IntKi), intent(out) :: ErrStat + character(*), intent(out) :: ErrMsg + + character(*), parameter :: RoutineName = 'amrex_find_subvols' character(c_char), allocatable :: dir_path(:) character(c_char), allocatable :: start_index(:) + integer(c_int), allocatable :: dir_indices_c(:) real(c_double) :: dt_c - integer(c_int) :: num_step, first_index, index_delta, subvol_c + integer(c_int) :: num_step, subvol_c integer(c_int) :: err_stat_c character(c_char) :: err_msg_c(ErrMsgLen) - integer(IntKi) :: i + integer(IntKi) :: i, stat + + ErrStat = ErrID_None + ErrMsg = "" #ifdef ENABLE_AMREX_LIB + if (NumStep < 1) then + call SetErrStat(ErrID_Fatal, "number of time steps must be at least 1", ErrStat, ErrMsg, RoutineName) + return + end if + ! Convert directory path to C type allocate(dir_path(len_trim(DirPath) + 1)) dir_path = transfer(trim(DirPath) // c_null_char, dir_path) @@ -181,20 +222,71 @@ subroutine amrex_find_subvols(DirPath, SubVol, DT, NumStep, StartIndex, & num_step = int(NumStep, c_int) dt_c = real(DT, c_double) - ! Call C++ routine to find first index and index delta values - call amrex_find_subvols_c(dir_path, subvol_c, dt_c, num_step, start_index, first_index, & - index_delta, err_stat_c, err_msg_c, ErrMsgLen) + ! Receive buffer for the C call. Allocated 1-based here and rebased to 0 on the way out, + ! since callers index the result with the (0-based) FAST.Farm time step number. + allocate(dir_indices_c(NumStep), stat=stat) + if (stat /= 0) then + call SetErrStat(ErrID_Fatal, "error allocating dir_indices_c", ErrStat, ErrMsg, RoutineName) + return + end if + dir_indices_c = -1_c_int + + ! Call C++ routine to build the time step -> directory index table + call amrex_find_subvols_c(dir_path, subvol_c, dt_c, num_step, start_index, dir_indices_c, & + err_stat_c, err_msg_c, ErrMsgLen) ! Transfer outputs to fortran types - FirstIndex = int(first_index, IntKi) - IndexDelta = int(index_delta, IntKi) ErrStat = int(err_stat_c, IntKi) ErrMsg = transfer(err_msg_c, ErrMsg) i = index(ErrMsg, c_null_char) if (i > 0) ErrMsg = ErrMsg(:i-1) + ! Leave DirIndices unallocated on failure so a caller that ignores ErrStat trips the + ! allocated() guard in ReadWindAMReX rather than silently reading zeros. + if (ErrStat >= AbortErrLev) return + + allocate(DirIndices(0:NumStep-1), stat=stat) + if (stat /= 0) then + call SetErrStat(ErrID_Fatal, "error allocating DirIndices", ErrStat, ErrMsg, RoutineName) + return + end if + DirIndices = int(dir_indices_c, IntKi) + #else - call SetErrStat(ErrID_Fatal, "AMReX library unavailable. Enable with -DAMREX_READER during compile with cmake on Linux, or change FAST.Farm Mod_AmbWind type", ErrStat, ErrMsg, "amrex_find_subvols") + call SetErrStat(ErrID_Fatal, "AMReX library unavailable. Enable with -DAMREX_READER during compile with cmake on Linux, or change FAST.Farm Mod_AmbWind type", ErrStat, ErrMsg, RoutineName) +#endif +end subroutine + +! Parse a plotfile Header directly, without AMReX, and return the grid information it holds. +! This is the fast path the sub-volume search uses once it has verified, on the starting +! directory, that the text agrees with amrex_read_header. Exposed so that agreement can be +! tested on real plotfiles. Ok is .false. if the header could not be parsed or describes a +! plotfile this reader does not accept. +subroutine amrex_parse_header_text(DirPath, Ok, Time, nXYZ, dXYZ, oXYZ) + character(*), intent(in) :: DirPath + logical, intent(out) :: Ok + real(DbKi), intent(out) :: Time + integer(IntKi), intent(out) :: nXYZ(3) + real(ReKi), intent(out) :: dXYZ(3) + real(ReKi), intent(out) :: oXYZ(3) + + character(c_char), allocatable :: dir_path(:) + integer(c_int) :: ok_c, dims(3) + real(c_double) :: t, origin(3), gridSpacing(3) + + Ok = .false.; Time = 0.0_DbKi; nXYZ = 0; dXYZ = 0.0_ReKi; oXYZ = 0.0_ReKi + +#ifdef ENABLE_AMREX_LIB + allocate(dir_path(len_trim(DirPath) + 1)) + dir_path = transfer(trim(DirPath) // c_null_char, dir_path) + call amrex_header_text_c(dir_path, t, dims, gridSpacing, origin, ok_c) + Ok = (ok_c /= 0) + if (Ok) then + Time = real(t, DbKi) + nXYZ = int(dims, IntKi) + dXYZ = real(gridSpacing, ReKi) + oXYZ = real(origin, ReKi) + end if #endif end subroutine diff --git a/modules/awae/src/amrex_utils.cpp b/modules/awae/src/amrex_utils.cpp index 0725ed0a91..9ac3f11d3a 100644 --- a/modules/awae/src/amrex_utils.cpp +++ b/modules/awae/src/amrex_utils.cpp @@ -6,6 +6,8 @@ #include #include #include +#include +#include #include @@ -69,8 +71,195 @@ void get_grid_bounds(const PlotFileData &pf, int level, // Define the variable names // const std::array var_names{"x_velocity", "y_velocity", "z_velocity"}; +// Extract the trailing directory index from a sub-volume path, e.g. "ffboxes_1_031150" -> 31150. +// Returns false if the suffix after the final '_' is not a plain non-negative integer, so that +// unrelated sibling directories (backups, renamed copies) are skipped rather than aborting the run. +// NOTE: the index must be compared numerically, never lexicographically: AMReX pads to a *minimum* +// width, so once a run passes 99999 steps a six-digit index sorts before a five-digit one as text. +bool parse_dir_index(const std::string &path, long long &index) +{ + const auto pos = path.find_last_of('_'); + if (pos == std::string::npos) + { + return false; + } + + const auto suffix = path.substr(pos + 1); + if (suffix.empty() || + suffix.find_first_not_of("0123456789") != std::string::npos) + { + return false; + } + + try + { + index = std::stoll(suffix); + } + catch (...) + { + return false; + } + + return true; +} + +// Grid metadata for one plotfile, as needed by the sub-volume search. +struct HeaderInfo +{ + double time{0.0}; + std::array dims{}; + std::array dx{}; + std::array origin{}; + int level_steps{-1}; +}; + +// Read a single-level plotfile Header directly, without constructing a PlotFileData. +// +// PlotFileData additionally opens Level_0/Cell_H and builds a DistributionMapping, which is +// wasted work when all that is wanted is the grid metadata; the sub-volume search does this once +// per directory, so on a large dataset the difference is hours. Everything needed is in the +// Header text: +// +// 1 version +// 2 ncomp +// 3 .. 2+ncomp variable names +// 3+ncomp spacedim +// 4+ncomp time +// 5+ncomp finest_level +// 6+ncomp prob_lo +// 7+ncomp prob_hi +// 8+ncomp ref_ratio (blank when finest_level is 0) +// 9+ncomp domain box, "((lo) (hi) (typ))" +// 10+ncomp level_steps (equals the directory index suffix) +// 11+ncomp cell size +// +// Returns false if anything does not parse -- or if the plotfile is not one this reader accepts +// (exactly three components, single level, level_steps matching the directory index when one is +// given) -- so the caller falls back to amrex_read_header_c, which reports the problem properly +// instead of guessing. Touches no AMReX global state. +bool parse_header_text(const std::string &dir, HeaderInfo &info, long long expect_index = -1) +{ + std::ifstream f(dir + "/Header"); + if (!f) + { + return false; + } + + std::vector line; + std::string s; + while (std::getline(f, s)) + { + line.push_back(s); + if (line.size() > 64) // everything of interest is near the top + { + break; + } + } + + try + { + if (line.size() < 2 || line[0].rfind("HyperCLaw", 0) != 0) + { + return false; + } + + // The reader takes the first three components positionally as velocity and requires + // exactly three. Anything else must be rejected loudly by amrex_read_header_c, not read. + const int ncomp = std::stoi(line[1]); + if (ncomp != 3) + { + return false; + } + + // 1-based line number -> 0-based index + const auto at = [&](int n) -> const std::string & { return line.at(n - 1); }; + + if (std::stoi(at(3 + ncomp)) != 3) // spacedim + { + return false; + } + info.time = std::stod(at(4 + ncomp)); + if (std::stoi(at(5 + ncomp)) != 0) // finest_level; the reader requires single level + { + return false; + } + + { + std::istringstream is(at(6 + ncomp)); + if (!(is >> info.origin[0] >> info.origin[1] >> info.origin[2])) + { + return false; + } + } + + // Domain box: "((0,0,0) (527,471,30) (0,0,0))" -> lo and hi index triples + { + auto b = at(9 + ncomp); + std::replace_if(b.begin(), b.end(), [](char c) { return c == '(' || c == ')' || c == ','; }, ' '); + std::istringstream is(b); + std::array lo{}, hi{}; + if (!(is >> lo[0] >> lo[1] >> lo[2] >> hi[0] >> hi[1] >> hi[2])) + { + return false; + } + + // level_steps is the solver step at which the plotfile was written and is what the + // directory suffix is generated from; a mismatch means the directory name does not + // describe its contents, so let the authoritative reader deal with it. + info.level_steps = std::stoi(at(10 + ncomp)); + if (expect_index >= 0 && info.level_steps != expect_index) + { + return false; + } + + std::istringstream ds(at(11 + ncomp)); + if (!(ds >> info.dx[0] >> info.dx[1] >> info.dx[2])) + { + return false; + } + + for (auto i = 0; i < 3; ++i) + { + if (hi[i] < lo[i]) + { + return false; + } + info.dims[i] = hi[i] - lo[i] + 1; + // Match amrex_read_header_c: problem origin + (grid index + 1/2) * cell size + info.origin[i] += (static_cast(lo[i]) + 0.5) * info.dx[i]; + } + } + } + catch (...) + { + return false; + } + + return true; +} + extern "C" { + // Parse the plotfile Header text directly (the fast path used by the sub-volume search) and + // return the grid information it yields. `ok` is 1 if the parse succeeded, 0 otherwise. This + // exists so the text parser can be tested against amrex_read_header_c on real plotfiles. + void amrex_header_text_c(char const *dir, double &time, int dims[3], double dx[3], double origin[3], int &ok) + { + HeaderInfo info; + ok = parse_header_text(std::string{dir}, info) ? 1 : 0; + if (ok == 0) + { + return; + } + time = info.time; + for (auto i = 0; i < 3; ++i) + { + dims[i] = info.dims[i]; + dx[i] = info.dx[i]; + origin[i] = info.origin[i]; + } + } + // Read the header information for the AMReX grid and return it. void amrex_read_header_c(char const *dir, double &time, int dims[3], double dx[3], double origin[3], int &err_stat, char *err_msg, int &err_msg_len) @@ -150,7 +339,13 @@ extern "C" // Read the XYZ velocity grid data into the FAST.Farm ambient wind data array [XYZ,NX,NY,NZ]. // This function cannot be called in parallel due to internal restrictions of the AMReX library. - void amrex_read_data_c(char const *dir, float *data, int &err_stat, char *err_msg, int &err_msg_len) + // + // `dims_expected` is the [NX,NY,NZ] extent of the caller's array. The plotfile must describe + // exactly that grid and must cover every cell of it: `data` is written by grid index, so a + // larger plotfile would write past the end of the caller's array, and one whose boxes leave + // holes would leave part of it holding whatever it held before. + void amrex_read_data_c(char const *dir, float *data, int const dims_expected[3], + int &err_stat, char *err_msg, int &err_msg_len) { const std::string routine{"amrex_read_data_c"}; @@ -180,6 +375,15 @@ extern "C" return; } + // Exactly three components are read positionally below; fewer would index past the FAB + const int ncomp = pf->nComp(); + if (ncomp != 3) + { + set_err(ErrID_Fatal, std::string{dir} + ": data dimensionality must be 3, got " + std::to_string(ncomp), + routine, err_stat, err_msg, err_msg_len); + return; + } + // Get overall grid bounds std::array dims{0}, gridLo{0}, gridHi{0}; get_grid_bounds(*pf, fine_level, gridLo, gridHi); @@ -188,15 +392,50 @@ extern "C" dims[i] = gridHi[i] - gridLo[i] + 1; } - // Get the variable names - const auto var_names = pf->varNames(); + // The grid must be exactly the one the caller allocated for. amrex_find_subvols_c checks + // this at initialization for every directory it matches, but it may have used the Header + // text fast path, which reads the domain box rather than the box array; re-check here + // against the destination array so a plotfile that disagrees can never be written out of + // bounds or leave the destination partly unwritten. + if ((dims[0] != dims_expected[0]) || (dims[1] != dims_expected[1]) || (dims[2] != dims_expected[2])) + { + const auto dims_str = "(" + std::to_string(dims[0]) + ", " + std::to_string(dims[1]) + ", " + std::to_string(dims[2]) + ")"; + const auto want_str = "(" + std::to_string(dims_expected[0]) + ", " + std::to_string(dims_expected[1]) + ", " + std::to_string(dims_expected[2]) + ")"; + set_err(ErrID_Fatal, std::string{dir} + ": grid dimensions " + dims_str + " do not match the " + want_str + + " grid established during initialization", + routine, err_stat, err_msg, err_msg_len); + return; + } - // Loop through variables - for (int ivar = 0; ivar < 3; ++ivar) + // The boxes of a plotfile box array are disjoint, so they tile the grid exactly when their + // volumes sum to its volume. Anything less leaves cells of `data` unwritten. { - // Get data for variable at given level - const auto &mf = pf->get(fine_level, var_names[ivar]); + const auto ba = pf->boxArray(fine_level); + long long covered = 0; + for (auto i = 0; i < ba.size(); ++i) + { + covered += static_cast(ba[i].numPts()); + } + const auto total = static_cast(dims[0]) * dims[1] * dims[2]; + if (covered != total) + { + set_err(ErrID_Fatal, std::string{dir} + ": box array covers " + std::to_string(covered) + " of " + + std::to_string(total) + " grid cells; the plotfile does not tile its grid", + routine, err_stat, err_msg, err_msg_len); + return; + } + } + + // Read every component in one pass. The per-variable overload re-reads the level once per + // variable, and the cost of a read is dominated by the number of boxes rather than by the + // volume of data, so three passes cost three times as much. Measured on a low-resolution + // sub-volume written with 93456 boxes: 65 s for three named reads against 17 s for one. + const auto &mf = pf->get(fine_level); + // Components are taken positionally, matching amrex_read_header_c's requirement that the + // first three are the X, Y and Z velocity in that order. + for (int ivar = 0; ivar < 3; ++ivar) + { // Loop through boxes of data for (MFIter mfi(mf); mfi.isValid(); ++mfi) { @@ -225,7 +464,7 @@ extern "C" { const auto gi = i - gridLo[0]; const auto di = get_grid_data_index(ivar, gi, gj, gk, 3, dims[0], dims[1]); - const auto v = fab(i, j, k); + const auto v = fab(i, j, k, ivar); data[di] = static_cast(v); } } @@ -234,25 +473,42 @@ extern "C" } } - // Search for AMReX directories based on given directory prefix, subvolume number, time step, total number of steps, - // and the starting index string (e.g. `00000`). This function returns first index as a number, and the delta - // between successive directory indices. It also checks that a sufficient number of directories are available, - // that the data matches the requested time step, and the grid properties are consistent (size, origin, spacing). + // Search for AMReX plotfile directories matching the given prefix and sub-volume number, and + // return the directory index to use for each of the `num_steps` requested time steps. + // + // Directories are matched to time steps by the simulation time recorded in each plotfile + // Header, NOT by any assumed stride between directory indices: the step claimed by a + // directory is round((header_time - start_time)/dt). This supports precursor data written + // with a varying solver time step -- for example an AMR-Wind run that transitions from + // time.initial_dt to fixed_dt, where the index stride changes but the output interval in + // time does not. + // + // Every step in [0, num_steps) must be claimed by a directory; a step claimed by none + // (missing data) is a fatal error, as is a step claimed twice (e.g. overlapping output from a + // restart) among the directories actually scanned. Note that the duplicate check is best + // effort rather than exhaustive: the walk stops once the table is complete and a directory + // lands past the window, so a stale duplicate at a higher index than that one is never read. + // Grid properties (size, origin, spacing) must be consistent across all steps. + // + // `dir_indices` must point to storage for at least num_steps ints. void amrex_find_subvols_c(char const *dir_prefix, int &subvol, double &dt, int &num_steps, char const *start_index, - int &first_index, int &index_delta, int &err_stat, char *err_msg, int &err_msg_len) + int *dir_indices, int &err_stat, char *err_msg, int &err_msg_len) { const std::string routine{"amrex_find_subvols_c"}; // Initialize error status and message to no error set_err(ErrID_None, "", routine, err_stat, err_msg, err_msg_len); + if (num_steps < 1) + { + set_err(ErrID_Fatal, "number of time steps must be at least 1, got " + std::to_string(num_steps), + routine, err_stat, err_msg, err_msg_len); + return; + } + // Construct path prefix based on directory prefix and subvolume number const std::filesystem::path path_prefix{std::string{dir_prefix} + "_" + std::to_string(subvol) + "_"}; - // Vector of index strings that match directory prefix and are - // greater than or equal to starting index - std::vector indices; - //---------------------------------------------------------------------- // Starting subvolume path //---------------------------------------------------------------------- @@ -263,7 +519,7 @@ extern "C" // If file does not exist, return error if (!std::filesystem::exists(first_path)) { - set_err(ErrID_Fatal, first_path + ": directory does not exist", + set_err(ErrID_Fatal, std::filesystem::absolute(first_path).string() + ": directory does not exist", routine, err_stat, err_msg, err_msg_len); return; } @@ -279,23 +535,111 @@ extern "C" return; } - // Save integer value of start index - first_index = std::stoi(start_index); + // The remaining directories are read with the much cheaper text parse. Confirm on this one + // directory that it agrees with the authoritative reader before trusting it for the rest. + // The two can in principle disagree: amrex_read_header_c takes the union of the box array + // in Level_0/Cell_H, whereas the Header records the domain box. They coincide for a + // single-level plotfile whose boxes tile its geometry, which is what the sub-volume writer + // produces -- but if that ever stops holding, fail loudly here rather than silently + // mismatching every subsequent directory. + bool use_fast_header = false; + { + HeaderInfo probe; + if (parse_header_text(first_path, probe)) + { + use_fast_header = (probe.dims == start_dims) && + (std::abs(probe.time - start_time) <= 1e-9 * std::max(1.0, std::abs(start_time))); + for (auto i = 0; i < 3 && use_fast_header; ++i) + { + use_fast_header = (std::abs(probe.dx[i] - start_dx[i]) <= 1e-8) && + (std::abs(probe.origin[i] - start_origin[i]) <= 1e-8); + } + } + } - // Add first index to list of valid indices - indices.emplace_back(first_index); + // Save integer value of start index + long long first_index_num{0}; + if (!parse_dir_index(first_path, first_index_num)) + { + set_err(ErrID_Fatal, std::string{start_index} + ": starting directory index must be a non-negative integer", + routine, err_stat, err_msg, err_msg_len); + return; + } + if (first_index_num > std::numeric_limits::max()) + { + set_err(ErrID_Fatal, std::string{start_index} + ": directory index exceeds the 32-bit range of the index table", + routine, err_stat, err_msg, err_msg_len); + return; + } //---------------------------------------------------------------------- - // Subsequent subvolume paths + // Time step matching tolerance //---------------------------------------------------------------------- - // If path prefix has parent directory use it, otherwise assume current directory - const auto parent_path = path_prefix.has_parent_path() ? path_prefix.parent_path() : "."; + // Tolerance on how far a directory's header time may sit from an exact multiple of dt. + // The error being absorbed is the drift a solver accumulates by summing its time step, + // which grows with absolute simulated time -- a precursor restarted at t = 3e4 s carries + // far more of it than one starting at zero -- so the tolerance is relative, with an + // absolute floor that preserves the historical behavior for runs starting near t = 0. + const auto step_tol = [&](double t) { + return std::max(1.0e-6, 1.0e-9 * (std::abs(start_time) + std::abs(t))); + }; + + // If the tolerance is an appreciable fraction of dt, step assignment is ambiguous and the + // caller should be told rather than silently given a clamped tolerance. + if (step_tol(start_time + static_cast(num_steps) * dt) >= 0.25 * dt) + { + set_err(ErrID_Fatal, path_prefix.string() + ": time step (" + std::to_string(dt) + + " s) is too small relative to the simulation time (" + std::to_string(start_time) + + " s) to identify time steps unambiguously", + routine, err_stat, err_msg, err_msg_len); + return; + } + + //---------------------------------------------------------------------- + // Assign each directory to the time step its header time corresponds to + //---------------------------------------------------------------------- - // Calculate maximum length of time from start time - const auto max_time = static_cast(num_steps) * dt; + // Directory index claiming each step, -1 if unclaimed + std::vector idx_of_step(num_steps, -1); + std::vector path_of_step(num_steps); + std::vector time_of_step(num_steps, 0.0); - // Loop through entries in the parent directory + // Closest directory that failed the residual test for each step, kept for diagnostics: + // when a step ends up unclaimed this is usually the file the user expected to fill it. + struct NearMiss + { + bool have{false}; + std::string path; + double time{0.0}; + double resid{0.0}; + double tol{0.0}; + }; + std::vector near_miss(num_steps); + + int n_before_start{0}, n_beyond_window{0}; + + // Seed step 0 from the start directory + idx_of_step[0] = first_index_num; + path_of_step[0] = first_path; + time_of_step[0] = start_time; + + // If path prefix has parent directory use it, otherwise assume current directory. The + // directory is scanned by name, so keep the final component of the prefix separately: an + // iterator over "." yields "./name", which does not begin with a prefix that has no + // directory of its own. + const auto parent_path = path_prefix.has_parent_path() ? path_prefix.parent_path() : std::filesystem::path{"."}; + const auto prefix_name = path_prefix.filename().string(); + + // Collect the candidate directories in one pass, then walk them in ascending index order. + // Ordering matters for cost, not correctness: within one run simulation time rises with the + // step counter, so once every step is claimed and a directory lands past the window, the + // remaining directories cannot add anything and the walk stops. Without that the search + // reads a header for every directory the LES ever wrote, however short the FAST.Farm run. + // The stop is conditional on the table being complete: leftovers from an earlier run with + // a different time step can put a later time on a lower index, and stopping on one of + // those would skip valid data that sorts after it. + std::vector> candidates; for (auto const &dir_entry : std::filesystem::directory_iterator{parent_path}) { // If entry is not a directory, continue @@ -307,48 +651,110 @@ extern "C" // Convert entry to path string const auto dir_path{dir_entry.path().string()}; - // If path doesn't contain the prefix, continue - if (dir_path.find(path_prefix.string()) == std::string::npos) + // If the name doesn't start with the prefix, continue. Anchored at position 0 so that + // an unrelated directory merely containing the prefix is not picked up. Compared on the + // final component only: when the prefix carries no directory of its own the iterator + // still yields entries as "./name", which no anchored compare against a bare prefix + // would ever match. + if (dir_entry.path().filename().string().rfind(prefix_name, 0) != 0) { continue; } - // Get the index string - const auto index = dir_path.substr(dir_path.find_last_of("_") + 1); + // Get the index, skipping entries whose suffix is not a plain integer + long long index{0}; + if (!parse_dir_index(dir_path, index)) + { + continue; + } - // If index string is less than starting index, continue - if (index.compare(start_index) <= 0) + // If index is not greater than the starting index, continue. This comparison must be + // numeric: a lexicographic compare drops every index wider than the starting index + // (e.g. "100030" sorts before "27150"), silently discarding data that is present. + if (index <= first_index_num) { continue; } + candidates.emplace_back(index, dir_path); + } + + std::sort(candidates.begin(), candidates.end()); + + // Step 0 is already claimed by the start directory + int n_claimed = 1; + + std::size_t visited = 0; + for (auto const &cand : candidates) + { + ++visited; + const auto index = cand.first; + const auto &dir_path = cand.second; + // Read the header double time{0.}; std::array dims; std::array dx, origin; - amrex_read_header_c(dir_path.c_str(), time, dims.data(), - dx.data(), origin.data(), err_stat, err_msg, err_msg_len); - if (err_stat != ErrID_None) + HeaderInfo hdr; + if (use_fast_header && parse_header_text(dir_path, hdr, index)) { - return; + time = hdr.time; + dims = hdr.dims; + dx = hdr.dx; + origin = hdr.origin; + } + else + { + amrex_read_header_c(dir_path.c_str(), time, dims.data(), + dx.data(), origin.data(), err_stat, err_msg, err_msg_len); + if (err_stat != ErrID_None) + { + return; + } } - // Get time delta from start time - const auto delta_time{time - start_time}; + const auto delta_time = time - start_time; + const auto tol = step_tol(time); - // If the delta time is greater than the max time plus dt (for safety), continue - if (delta_time > (max_time + dt / 4.)) + if (delta_time < -tol) { + ++n_before_start; continue; } - // If delta time is not a multiple of dt, continue - // (remainder must be nearly zero or nearly equal to dt) - // A tolerance of 1e-6 seconds seems reasonable - const auto remainder{std::fmod(delta_time, dt)}; - if (!((std::abs(remainder - 0.0) <= 1e-6) || - ((std::abs(remainder - dt) <= 1e-6)))) + // Nearest time step, and how far this directory sits from it + const auto step = std::lround(delta_time / dt); + const auto resid = std::abs(delta_time - static_cast(step) * dt); + + if (step < 0) + { + ++n_before_start; + continue; + } + if (step >= static_cast(num_steps)) { + ++n_beyond_window; + if (n_claimed == num_steps) + { + // Table complete and this directory is past the window: nothing after it in + // ascending index order can be needed. Count the rest as skipped and stop. + n_beyond_window += static_cast(candidates.size() - visited); + break; + } + // Something is still missing, so do not trust index order to imply time order; + // keep walking and let a genuinely missing step be reported after the full scan. + continue; + } + + // Not on a step boundary. This is how deliberately decimated output is skipped: when + // dt is a multiple of the file cadence, the intermediate files land here. + if (resid > tol) + { + auto &nm = near_miss[step]; + if (!nm.have || resid < nm.resid) + { + nm = NearMiss{true, dir_path, time, resid, tol}; + } continue; } @@ -379,64 +785,75 @@ extern "C" return; } - // Add index to list of indices - indices.emplace_back(std::stoi(index)); + // Two directories cannot represent the same instant in time + if (idx_of_step[step] >= 0) + { + std::string msg{path_prefix.string() + ": two sub-volume directories claim time step "}; + msg += std::to_string(step) + " (expected header time " + std::to_string(start_time + static_cast(step) * dt) + " s): '"; + msg += path_of_step[step] + "' (header t = " + std::to_string(time_of_step[step]) + " s) and '"; + msg += dir_path + "' (header t = " + std::to_string(time) + " s). Each time step must be represented by "; + msg += "exactly one directory; this usually means output from two different runs (e.g. a restart that "; + msg += "re-wrote overlapping times) is present. Remove or move the stale directories."; + set_err(ErrID_Fatal, msg, routine, err_stat, err_msg, err_msg_len); + return; + } + + idx_of_step[step] = index; + path_of_step[step] = dir_path; + time_of_step[step] = time; + ++n_claimed; } //---------------------------------------------------------------------- - // Check indices + // Every step must be accounted for //---------------------------------------------------------------------- - // Check that more than one index was found - if (indices.size() < 2) + for (int s = 0; s < num_steps; ++s) { - set_err(ErrID_Fatal, path_prefix.string() + ": only 1 subvolume found, at least 2 required", - routine, err_stat, err_msg, err_msg_len); - return; - } + if (idx_of_step[s] >= 0) + { + if (idx_of_step[s] > std::numeric_limits::max()) + { + set_err(ErrID_Fatal, path_of_step[s] + ": directory index exceeds the 32-bit range of the index table", + routine, err_stat, err_msg, err_msg_len); + return; + } + dir_indices[s] = static_cast(idx_of_step[s]); + continue; + } - // Sort indices in ascending order - std::sort(indices.begin(), indices.end()); + const auto want_time = start_time + static_cast(s) * dt; - // If more indices found that requested steps, discard extra - if (indices.size() > num_steps) - { - indices.resize(num_steps); - } - - // Calculate delta between first two indices - const auto first_index_delta = indices[1] - indices[0]; + std::string msg{path_prefix.string() + ": no sub-volume directory was found for time step "}; + msg += std::to_string(s) + " of " + std::to_string(num_steps) + ". Expected header time "; + msg += std::to_string(want_time) + " s = " + std::to_string(start_time) + " s (start directory '"; + msg += first_path + "') + " + std::to_string(s) + " * dt (" + std::to_string(dt) + " s), matched to "; + msg += "within " + std::to_string(step_tol(want_time)) + " s."; - // If fewer indices found than requested steps, return error - if (indices.size() < num_steps) - { - std::string msg{path_prefix.string() + ": "}; - msg += "found " + std::to_string(indices.size()) + " dirs, "; - msg += std::to_string(num_steps) + " dirs were requested "; - msg += "with a dt of " + std::to_string(dt) + " seconds ("; - msg += "first index=" + std::to_string(indices[0]) + ", "; - msg += "last index=" + std::to_string(indices.back()) + ", "; - msg += "index step=" + std::to_string(first_index_delta) + ")"; - set_err(ErrID_Fatal, msg, routine, err_stat, err_msg, err_msg_len); - return; - } + if (near_miss[s].have) + { + msg += " Directory '" + near_miss[s].path + "' exists with header time "; + msg += std::to_string(near_miss[s].time) + " s, which is " + std::to_string(near_miss[s].resid); + msg += " s (" + std::to_string(near_miss[s].resid / dt) + " * dt) from the expected time -- outside "; + msg += "the " + std::to_string(near_miss[s].tol) + " s tolerance applied to it."; + } - // Loop through indices and check that none are missing - // ie, same delta between all adjacent indicies - for (auto i = 1; i < indices.size(); ++i) - { - // Calculate index delta between current and previous indices - index_delta = indices[i] - indices[i - 1]; + msg += " Sub-volume directories are matched to FAST.Farm time steps by the simulation time recorded in "; + msg += "their Header; the directory index stride is irrelevant and may vary."; - // If delta doesn't match first delta, return error - if (index_delta != first_index_delta) + if (n_before_start > 0) { - std::string msg{path_prefix.string() + ": "}; - msg += "inconsistent delta between indices '" + std::to_string(indices[i - 1]); - msg += "' and '" + std::to_string(indices[i]) + "'"; - set_err(ErrID_Fatal, msg, routine, err_stat, err_msg, err_msg_len); - return; + msg += " (" + std::to_string(n_before_start) + " directories were skipped because their header time "; + msg += "precedes the start directory.)"; } + if (n_beyond_window > 0) + { + msg += " (" + std::to_string(n_beyond_window) + " directories were skipped because their header time "; + msg += "is past the end of the requested window.)"; + } + + set_err(ErrID_Fatal, msg, routine, err_stat, err_msg, err_msg_len); + return; } } } diff --git a/modules/awae/tests/data/subvolstale_0_00000/Header b/modules/awae/tests/data/subvolstale_0_00000/Header new file mode 100644 index 0000000000..033743c58c --- /dev/null +++ b/modules/awae/tests/data/subvolstale_0_00000/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +0.0 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +0 +1 1 1 +0 +0 +0 2 0.0 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolstale_0_00000/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolstale_0_00000/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolstale_0_00000/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolstale_0_00000/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolstale_0_00000/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolstale_0_00000/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolstale_0_00000/Level_0/Cell_H b/modules/awae/tests/data/subvolstale_0_00000/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolstale_0_00000/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolstale_0_00002/Header b/modules/awae/tests/data/subvolstale_0_00002/Header new file mode 100644 index 0000000000..3fc2e7d63a --- /dev/null +++ b/modules/awae/tests/data/subvolstale_0_00002/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +5.0 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +2 +1 1 1 +0 +0 +0 2 5.0 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolstale_0_00002/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolstale_0_00002/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolstale_0_00002/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolstale_0_00002/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolstale_0_00002/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolstale_0_00002/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolstale_0_00002/Level_0/Cell_H b/modules/awae/tests/data/subvolstale_0_00002/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolstale_0_00002/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolstale_0_00004/Header b/modules/awae/tests/data/subvolstale_0_00004/Header new file mode 100644 index 0000000000..b3f54c076d --- /dev/null +++ b/modules/awae/tests/data/subvolstale_0_00004/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +0.1 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +4 +1 1 1 +0 +0 +0 2 0.1 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolstale_0_00004/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolstale_0_00004/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolstale_0_00004/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolstale_0_00004/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolstale_0_00004/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolstale_0_00004/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolstale_0_00004/Level_0/Cell_H b/modules/awae/tests/data/subvolstale_0_00004/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolstale_0_00004/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolstale_0_00008/Header b/modules/awae/tests/data/subvolstale_0_00008/Header new file mode 100644 index 0000000000..155c842db8 --- /dev/null +++ b/modules/awae/tests/data/subvolstale_0_00008/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +0.2 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +8 +1 1 1 +0 +0 +0 2 0.2 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolstale_0_00008/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolstale_0_00008/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolstale_0_00008/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolstale_0_00008/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolstale_0_00008/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolstale_0_00008/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolstale_0_00008/Level_0/Cell_H b/modules/awae/tests/data/subvolstale_0_00008/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolstale_0_00008/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvoltiled_1_31220/Header b/modules/awae/tests/data/subvoltiled_1_31220/Header new file mode 100644 index 0000000000..1e30993c40 --- /dev/null +++ b/modules/awae/tests/data/subvoltiled_1_31220/Header @@ -0,0 +1,22 @@ +HyperCLaw-V1.1 +3 +velocity_maskedx +velocity_maskedy +velocity_maskedz +3 +29999.999999998108 +0 +681089.49589999998 5058948.6238000002 356.5 +681233.49589999998 5059092.6238000002 508.5 + +((0,0,0) (17,17,18) (0,0,0)) +31220 +8 8 8 +0 +0 +0 1 29999.999999998108 +31220 +681089.49589999998 681233.49589999998 +5058948.6238000002 5059092.6238000002 +356.5 508.5 +Level_0/Cell diff --git a/modules/awae/tests/data/subvoltiled_1_31220/Level_0/Cell_D_00018 b/modules/awae/tests/data/subvoltiled_1_31220/Level_0/Cell_D_00018 new file mode 100644 index 0000000000..5cb1c1f576 Binary files /dev/null and b/modules/awae/tests/data/subvoltiled_1_31220/Level_0/Cell_D_00018 differ diff --git a/modules/awae/tests/data/subvoltiled_1_31220/Level_0/Cell_H b/modules/awae/tests/data/subvoltiled_1_31220/Level_0/Cell_H new file mode 100644 index 0000000000..b6a7e16db4 --- /dev/null +++ b/modules/awae/tests/data/subvoltiled_1_31220/Level_0/Cell_H @@ -0,0 +1,16 @@ +1 +1 +3 +0 +(1 0 +((0,0,0) (17,17,18) (0,0,0)) +) +1 +FabOnDisk: Cell_D_00018 0 + +1,3 +1.75883225107680374e+00,-9.75127050900527426e-01,-1.54867083618545753e+00, + +1,3 +9.12948102152889085e+00,2.87040254526559790e+00,2.02065795062786080e+00, + diff --git a/modules/awae/tests/data/subvoltiled_1_31224/Header b/modules/awae/tests/data/subvoltiled_1_31224/Header new file mode 100644 index 0000000000..cda2397135 --- /dev/null +++ b/modules/awae/tests/data/subvoltiled_1_31224/Header @@ -0,0 +1,22 @@ +HyperCLaw-V1.1 +3 +velocity_maskedx +velocity_maskedy +velocity_maskedz +3 +30000.099999998114 +0 +681089.49589999998 5058948.6238000002 356.5 +681233.49589999998 5059092.6238000002 508.5 + +((0,0,0) (17,17,18) (0,0,0)) +31224 +8 8 8 +0 +0 +0 1 30000.099999998114 +31224 +681089.49589999998 681233.49589999998 +5058948.6238000002 5059092.6238000002 +356.5 508.5 +Level_0/Cell diff --git a/modules/awae/tests/data/subvoltiled_1_31224/Level_0/Cell_D_00018 b/modules/awae/tests/data/subvoltiled_1_31224/Level_0/Cell_D_00018 new file mode 100644 index 0000000000..5cb1c1f576 Binary files /dev/null and b/modules/awae/tests/data/subvoltiled_1_31224/Level_0/Cell_D_00018 differ diff --git a/modules/awae/tests/data/subvoltiled_1_31224/Level_0/Cell_H b/modules/awae/tests/data/subvoltiled_1_31224/Level_0/Cell_H new file mode 100644 index 0000000000..6b99a044d6 --- /dev/null +++ b/modules/awae/tests/data/subvoltiled_1_31224/Level_0/Cell_H @@ -0,0 +1,16 @@ +1 +1 +3 +0 +(1 0 +((0,0,0) (17,17,18) (0,0,0)) +) +1 +FabOnDisk: Cell_D_00018 0 + +1,3 +1.79769313486231571e+308,1.79769313486231571e+308,1.79769313486231571e+308, + +1,3 +-1.79769313486231571e+308,-1.79769313486231571e+308,-1.79769313486231571e+308, + diff --git a/modules/awae/tests/data/subvoltiled_1_31228/Header b/modules/awae/tests/data/subvoltiled_1_31228/Header new file mode 100644 index 0000000000..5fe57d6a2f --- /dev/null +++ b/modules/awae/tests/data/subvoltiled_1_31228/Header @@ -0,0 +1,22 @@ +HyperCLaw-V1.1 +3 +velocity_maskedx +velocity_maskedy +velocity_maskedz +3 +30000.19999999812 +0 +681089.49589999998 5058948.6238000002 356.5 +681233.49589999998 5059092.6238000002 508.5 + +((0,0,0) (17,17,18) (0,0,0)) +31228 +8 8 8 +0 +0 +0 1 30000.19999999812 +31228 +681089.49589999998 681233.49589999998 +5058948.6238000002 5059092.6238000002 +356.5 508.5 +Level_0/Cell diff --git a/modules/awae/tests/data/subvoltiled_1_31228/Level_0/Cell_D_00018 b/modules/awae/tests/data/subvoltiled_1_31228/Level_0/Cell_D_00018 new file mode 100644 index 0000000000..5cb1c1f576 Binary files /dev/null and b/modules/awae/tests/data/subvoltiled_1_31228/Level_0/Cell_D_00018 differ diff --git a/modules/awae/tests/data/subvoltiled_1_31228/Level_0/Cell_H b/modules/awae/tests/data/subvoltiled_1_31228/Level_0/Cell_H new file mode 100644 index 0000000000..6b99a044d6 --- /dev/null +++ b/modules/awae/tests/data/subvoltiled_1_31228/Level_0/Cell_H @@ -0,0 +1,16 @@ +1 +1 +3 +0 +(1 0 +((0,0,0) (17,17,18) (0,0,0)) +) +1 +FabOnDisk: Cell_D_00018 0 + +1,3 +1.79769313486231571e+308,1.79769313486231571e+308,1.79769313486231571e+308, + +1,3 +-1.79769313486231571e+308,-1.79769313486231571e+308,-1.79769313486231571e+308, + diff --git a/modules/awae/tests/data/subvolvardt_0_00000/Header b/modules/awae/tests/data/subvolvardt_0_00000/Header new file mode 100644 index 0000000000..a4be1c0417 --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_0_00000/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +29999.999999998545 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +0 +1 1 1 +0 +0 +0 2 29999.999999998545 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolvardt_0_00000/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolvardt_0_00000/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_0_00000/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolvardt_0_00000/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolvardt_0_00000/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_0_00000/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolvardt_0_00000/Level_0/Cell_H b/modules/awae/tests/data/subvolvardt_0_00000/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_0_00000/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolvardt_0_00004/Header b/modules/awae/tests/data/subvolvardt_0_00004/Header new file mode 100644 index 0000000000..1c48a394b0 --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_0_00004/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +30000.099999998543 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +4 +1 1 1 +0 +0 +0 2 30000.099999998543 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolvardt_0_00004/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolvardt_0_00004/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_0_00004/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolvardt_0_00004/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolvardt_0_00004/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_0_00004/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolvardt_0_00004/Level_0/Cell_H b/modules/awae/tests/data/subvolvardt_0_00004/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_0_00004/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolvardt_0_00008/Header b/modules/awae/tests/data/subvolvardt_0_00008/Header new file mode 100644 index 0000000000..c462b0d577 --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_0_00008/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +30000.199999998546 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +8 +1 1 1 +0 +0 +0 2 30000.199999998546 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolvardt_0_00008/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolvardt_0_00008/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_0_00008/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolvardt_0_00008/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolvardt_0_00008/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_0_00008/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolvardt_0_00008/Level_0/Cell_H b/modules/awae/tests/data/subvolvardt_0_00008/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_0_00008/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolvardt_0_00010/Header b/modules/awae/tests/data/subvolvardt_0_00010/Header new file mode 100644 index 0000000000..f40dcd9822 --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_0_00010/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +30000.300000004358 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +10 +1 1 1 +0 +0 +0 2 30000.300000004358 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolvardt_0_00010/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolvardt_0_00010/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_0_00010/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolvardt_0_00010/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolvardt_0_00010/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_0_00010/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolvardt_0_00010/Level_0/Cell_H b/modules/awae/tests/data/subvolvardt_0_00010/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_0_00010/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolvardt_0_00012/Header b/modules/awae/tests/data/subvolvardt_0_00012/Header new file mode 100644 index 0000000000..05a1b0998f --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_0_00012/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +30000.400000014357 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +12 +1 1 1 +0 +0 +0 2 30000.400000014357 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolvardt_0_00012/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolvardt_0_00012/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_0_00012/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolvardt_0_00012/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolvardt_0_00012/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_0_00012/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolvardt_0_00012/Level_0/Cell_H b/modules/awae/tests/data/subvolvardt_0_00012/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_0_00012/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolvardt_1_00000/Header b/modules/awae/tests/data/subvolvardt_1_00000/Header new file mode 100644 index 0000000000..b091acde08 --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_1_00000/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +29999.999999998545 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +0 +1 1 1 +0 +0 +0 2 29999.999999998545 +0 +6 8 +6 9 +6 9 +6 8 +6 9 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolvardt_1_00000/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolvardt_1_00000/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..2053e741d5 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_1_00000/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolvardt_1_00000/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolvardt_1_00000/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..1d62f2a852 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_1_00000/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolvardt_1_00000/Level_0/Cell_H b/modules/awae/tests/data/subvolvardt_1_00000/Level_0/Cell_H new file mode 100644 index 0000000000..6f35cc6918 --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_1_00000/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (7,8,8) (0,0,0)) +((6,6,9) (7,8,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolvardt_1_00004/Header b/modules/awae/tests/data/subvolvardt_1_00004/Header new file mode 100644 index 0000000000..c01b4c614a --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_1_00004/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +30000.099999998543 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +4 +1 1 1 +0 +0 +0 2 30000.099999998543 +0 +6 8 +6 9 +6 9 +6 8 +6 9 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolvardt_1_00004/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolvardt_1_00004/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..2053e741d5 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_1_00004/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolvardt_1_00004/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolvardt_1_00004/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..1d62f2a852 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_1_00004/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolvardt_1_00004/Level_0/Cell_H b/modules/awae/tests/data/subvolvardt_1_00004/Level_0/Cell_H new file mode 100644 index 0000000000..6f35cc6918 --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_1_00004/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (7,8,8) (0,0,0)) +((6,6,9) (7,8,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolvardt_1_00005/Header b/modules/awae/tests/data/subvolvardt_1_00005/Header new file mode 100644 index 0000000000..6f801a696f --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_1_00005/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +30000.099999998543 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +5 +1 1 1 +0 +0 +0 2 30000.099999998543 +0 +6 8 +6 9 +6 9 +6 8 +6 9 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolvardt_1_00005/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolvardt_1_00005/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..2053e741d5 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_1_00005/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolvardt_1_00005/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolvardt_1_00005/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..1d62f2a852 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_1_00005/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolvardt_1_00005/Level_0/Cell_H b/modules/awae/tests/data/subvolvardt_1_00005/Level_0/Cell_H new file mode 100644 index 0000000000..6f35cc6918 --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_1_00005/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (7,8,8) (0,0,0)) +((6,6,9) (7,8,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolvardt_2_00000/Header b/modules/awae/tests/data/subvolvardt_2_00000/Header new file mode 100644 index 0000000000..a4be1c0417 --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_2_00000/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +29999.999999998545 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +0 +1 1 1 +0 +0 +0 2 29999.999999998545 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolvardt_2_00000/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolvardt_2_00000/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_2_00000/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolvardt_2_00000/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolvardt_2_00000/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_2_00000/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolvardt_2_00000/Level_0/Cell_H b/modules/awae/tests/data/subvolvardt_2_00000/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_2_00000/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolvardt_2_00004/Header b/modules/awae/tests/data/subvolvardt_2_00004/Header new file mode 100644 index 0000000000..e9e30eccb3 --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_2_00004/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +30000.101 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +4 +1 1 1 +0 +0 +0 2 30000.101 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolvardt_2_00004/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolvardt_2_00004/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_2_00004/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolvardt_2_00004/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolvardt_2_00004/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolvardt_2_00004/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolvardt_2_00004/Level_0/Cell_H b/modules/awae/tests/data/subvolvardt_2_00004/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolvardt_2_00004/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolwide_0_100000/Header b/modules/awae/tests/data/subvolwide_0_100000/Header new file mode 100644 index 0000000000..d81671cacc --- /dev/null +++ b/modules/awae/tests/data/subvolwide_0_100000/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +0.1 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +100000 +1 1 1 +0 +0 +0 2 0.1 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolwide_0_100000/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolwide_0_100000/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolwide_0_100000/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolwide_0_100000/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolwide_0_100000/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolwide_0_100000/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolwide_0_100000/Level_0/Cell_H b/modules/awae/tests/data/subvolwide_0_100000/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolwide_0_100000/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolwide_0_100002/Header b/modules/awae/tests/data/subvolwide_0_100002/Header new file mode 100644 index 0000000000..b1160d7915 --- /dev/null +++ b/modules/awae/tests/data/subvolwide_0_100002/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +0.2 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +100002 +1 1 1 +0 +0 +0 2 0.2 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolwide_0_100002/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolwide_0_100002/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolwide_0_100002/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolwide_0_100002/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolwide_0_100002/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolwide_0_100002/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolwide_0_100002/Level_0/Cell_H b/modules/awae/tests/data/subvolwide_0_100002/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolwide_0_100002/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/data/subvolwide_0_99998/Header b/modules/awae/tests/data/subvolwide_0_99998/Header new file mode 100644 index 0000000000..23d680728f --- /dev/null +++ b/modules/awae/tests/data/subvolwide_0_99998/Header @@ -0,0 +1,25 @@ +HyperCLaw-V1.1 +3 +x_velocity +y_velocity +z_velocity +3 +0.0 +0 +0 0 0 +256 256 128 + +((0,0,0) (255,255,127) (0,0,0)) +99998 +1 1 1 +0 +0 +0 2 0.0 +0 +6 9 +6 10 +6 9 +6 9 +6 10 +9 11 +Level_0/Cell diff --git a/modules/awae/tests/data/subvolwide_0_99998/Level_0/Cell_D_00064 b/modules/awae/tests/data/subvolwide_0_99998/Level_0/Cell_D_00064 new file mode 100644 index 0000000000..fd33ae991c Binary files /dev/null and b/modules/awae/tests/data/subvolwide_0_99998/Level_0/Cell_D_00064 differ diff --git a/modules/awae/tests/data/subvolwide_0_99998/Level_0/Cell_D_00065 b/modules/awae/tests/data/subvolwide_0_99998/Level_0/Cell_D_00065 new file mode 100644 index 0000000000..93c211c471 Binary files /dev/null and b/modules/awae/tests/data/subvolwide_0_99998/Level_0/Cell_D_00065 differ diff --git a/modules/awae/tests/data/subvolwide_0_99998/Level_0/Cell_H b/modules/awae/tests/data/subvolwide_0_99998/Level_0/Cell_H new file mode 100644 index 0000000000..523080d05c --- /dev/null +++ b/modules/awae/tests/data/subvolwide_0_99998/Level_0/Cell_H @@ -0,0 +1,20 @@ +1 +1 +3 +0 +(2 0 +((6,6,6) (8,9,8) (0,0,0)) +((6,6,9) (8,9,10) (0,0,0)) +) +2 +FabOnDisk: Cell_D_00064 0 +FabOnDisk: Cell_D_00065 0 + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + +2,3 +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, +1.00000000000000000e+01,0.00000000000000000e+00,0.00000000000000000e+00, + diff --git a/modules/awae/tests/test_AMReX_reader.f90 b/modules/awae/tests/test_AMReX_reader.f90 index b382ed4750..b3833f8492 100644 --- a/modules/awae/tests/test_AMReX_reader.f90 +++ b/modules/awae/tests/test_AMReX_reader.f90 @@ -2,6 +2,8 @@ module test_AMReX_reader use testdrive, only: new_unittest, unittest_type, error_type, check use NWTC_Library use amrex_utils + use AWAE_Types, only: AWAE_ParameterType + use AWAE_IO, only: ReadWindAMReX implicit none @@ -14,7 +16,23 @@ subroutine run_test_AMReX_reader(testsuite) new_unittest("AMReX_test_read_subvol_1", AMReX_test_read_subvol_1), & new_unittest("AMReX_test_amrex_find_subvols_1", AMReX_test_amrex_find_subvols_1), & new_unittest("AMReX_test_amrex_find_subvols_2", AMReX_test_amrex_find_subvols_2), & - new_unittest("AMReX_test_amrex_find_subvols_3", AMReX_test_amrex_find_subvols_3) & + new_unittest("AMReX_test_amrex_find_subvols_3", AMReX_test_amrex_find_subvols_3), & + new_unittest("AMReX_test_find_subvols_variable_stride", AMReX_test_find_subvols_variable_stride), & + new_unittest("AMReX_test_find_subvols_missing_step", AMReX_test_find_subvols_missing_step), & + new_unittest("AMReX_test_find_subvols_duplicate_step", AMReX_test_find_subvols_duplicate_step), & + new_unittest("AMReX_test_find_subvols_out_of_tolerance", AMReX_test_find_subvols_out_of_tolerance), & + new_unittest("AMReX_test_find_subvols_wide_index", AMReX_test_find_subvols_wide_index), & + new_unittest("AMReX_test_ReadWindAMReX_lookup", AMReX_test_ReadWindAMReX_lookup), & + new_unittest("AMReX_test_ReadWindAMReX_out_of_range", AMReX_test_ReadWindAMReX_out_of_range), & + new_unittest("AMReX_test_ReadWindAMReX_wide_index", AMReX_test_ReadWindAMReX_wide_index), & + new_unittest("AMReX_test_ReadWindAMReX_no_table", AMReX_test_ReadWindAMReX_no_table), & + new_unittest("AMReX_test_find_subvols_stale_beyond_window", AMReX_test_find_subvols_stale_beyond_window), & + new_unittest("AMReX_test_find_subvols_real_tiled", AMReX_test_find_subvols_real_tiled), & + new_unittest("AMReX_test_read_real_tiled", AMReX_test_read_real_tiled), & + new_unittest("AMReX_test_header_text_agrees_on_tiled", AMReX_test_header_text_agrees_on_tiled), & + new_unittest("AMReX_test_header_text_disagrees_on_subset", AMReX_test_header_text_disagrees_on_subset), & + new_unittest("AMReX_test_find_subvols_prefix_form", AMReX_test_find_subvols_prefix_form), & + new_unittest("AMReX_test_read_data_wrong_dims", AMReX_test_read_data_wrong_dims) & ] end subroutine @@ -27,7 +45,8 @@ subroutine AMReX_test_read_subvol_0(error) character(ErrMsgLen) :: ErrMsg integer(IntKi) :: i, j, k integer(IntKi) :: dims(3) - real(ReKi) :: time, origin(3), gridSpacing(3), gridSize(3), bounds(2,3) + real(DbKi) :: time + real(ReKi) :: origin(3), gridSpacing(3), gridSize(3), bounds(2,3) ! Read header call amrex_read_header(trim(DirPath), time, dims, gridSpacing, origin, ErrStat, ErrMsg) @@ -35,14 +54,14 @@ subroutine AMReX_test_read_subvol_0(error) if (ErrStat /= ErrID_None) print*, "ErrMsg = ", ErrMsg ! Calculate grid size - gridSize = gridSpacing*real(dims - 1, c_double) + gridSize = gridSpacing*real(dims - 1, ReKi) ! Calculate the grid bounds bounds(1,:) = origin bounds(2,:) = origin + gridSize ! Check time - call check(error, time, 0.6_c_double); if (allocated(error)) return + call check(error, time, 0.6_DbKi); if (allocated(error)) return ! Check dimensions call check(error, dims(1), 3_c_int, more="dims(1)"); if (allocated(error)) return @@ -50,24 +69,24 @@ subroutine AMReX_test_read_subvol_0(error) call check(error, dims(3), 5_c_int, more="dims(3)"); if (allocated(error)) return ! Check spacing - call check(error, gridSpacing(1), 1.0_c_double, more="gridSpacing(1)"); if (allocated(error)) return - call check(error, gridSpacing(2), 1.0_c_double, more="gridSpacing(2)"); if (allocated(error)) return - call check(error, gridSpacing(3), 1.0_c_double, more="gridSpacing(3)"); if (allocated(error)) return + call check(error, gridSpacing(1), 1.0_ReKi, more="gridSpacing(1)"); if (allocated(error)) return + call check(error, gridSpacing(2), 1.0_ReKi, more="gridSpacing(2)"); if (allocated(error)) return + call check(error, gridSpacing(3), 1.0_ReKi, more="gridSpacing(3)"); if (allocated(error)) return ! Check grid size - call check(error, gridSize(1), 2.0_c_double, more="gridSize(1)"); if (allocated(error)) return - call check(error, gridSize(2), 3.0_c_double, more="gridSize(2)"); if (allocated(error)) return - call check(error, gridSize(3), 4.0_c_double, more="gridSize(3)"); if (allocated(error)) return + call check(error, gridSize(1), 2.0_ReKi, more="gridSize(1)"); if (allocated(error)) return + call check(error, gridSize(2), 3.0_ReKi, more="gridSize(2)"); if (allocated(error)) return + call check(error, gridSize(3), 4.0_ReKi, more="gridSize(3)"); if (allocated(error)) return ! Check lower bounds - call check(error, bounds(1,1), 6.5_c_double, more="origin(1)"); if (allocated(error)) return - call check(error, bounds(1,2), 6.5_c_double, more="origin(2)"); if (allocated(error)) return - call check(error, bounds(1,3), 6.5_c_double, more="origin(3)"); if (allocated(error)) return + call check(error, bounds(1,1), 6.5_ReKi, more="origin(1)"); if (allocated(error)) return + call check(error, bounds(1,2), 6.5_ReKi, more="origin(2)"); if (allocated(error)) return + call check(error, bounds(1,3), 6.5_ReKi, more="origin(3)"); if (allocated(error)) return ! Check upper bounds - call check(error, bounds(2,1), 8.5_c_double, more="ub(1)"); if (allocated(error)) return - call check(error, bounds(2,2), 9.5_c_double, more="ub(2)"); if (allocated(error)) return - call check(error, bounds(2,3), 10.5_c_double, more="ub(3)"); if (allocated(error)) return + call check(error, bounds(2,1), 8.5_ReKi, more="ub(1)"); if (allocated(error)) return + call check(error, bounds(2,2), 9.5_ReKi, more="ub(2)"); if (allocated(error)) return + call check(error, bounds(2,3), 10.5_ReKi, more="ub(3)"); if (allocated(error)) return ! Display grid properties print*, "dir = ", trim(DirPath) @@ -109,7 +128,8 @@ subroutine AMReX_test_read_subvol_1(error) character(ErrMsgLen) :: ErrMsg integer(IntKi) :: i, j, k integer(IntKi) :: dims(3) - real(ReKi) :: time, origin(3), gridSpacing(3), gridSize(3), bounds(2,3) + real(DbKi) :: time + real(ReKi) :: origin(3), gridSpacing(3), gridSize(3), bounds(2,3) ! Read header call amrex_read_header(trim(DirPath), time, dims, gridSpacing, origin, ErrStat, ErrMsg) @@ -117,14 +137,14 @@ subroutine AMReX_test_read_subvol_1(error) if (ErrStat /= ErrID_None) print*, "ErrMsg = ", ErrMsg ! Calculate grid size - gridSize = gridSpacing*real(dims - 1, c_double) + gridSize = gridSpacing*real(dims - 1, ReKi) ! Calculate the grid bounds bounds(1,:) = origin bounds(2,:) = origin + gridSize ! Check time - call check(error, time, 1.6_c_double); if (allocated(error)) return + call check(error, time, 1.6_DbKi); if (allocated(error)) return ! Check dimensions call check(error, dims(1), 2_c_int, more="dims(1)"); if (allocated(error)) return @@ -132,24 +152,24 @@ subroutine AMReX_test_read_subvol_1(error) call check(error, dims(3), 5_c_int, more="dims(3)"); if (allocated(error)) return ! Check spacing - call check(error, gridSpacing(1), 1.0_c_double, more="gridSpacing(1)"); if (allocated(error)) return - call check(error, gridSpacing(2), 1.0_c_double, more="gridSpacing(2)"); if (allocated(error)) return - call check(error, gridSpacing(3), 1.0_c_double, more="gridSpacing(3)"); if (allocated(error)) return + call check(error, gridSpacing(1), 1.0_ReKi, more="gridSpacing(1)"); if (allocated(error)) return + call check(error, gridSpacing(2), 1.0_ReKi, more="gridSpacing(2)"); if (allocated(error)) return + call check(error, gridSpacing(3), 1.0_ReKi, more="gridSpacing(3)"); if (allocated(error)) return ! Check grid size - call check(error, gridSize(1), 1.0_c_double, more="gridSize(1)"); if (allocated(error)) return - call check(error, gridSize(2), 2.0_c_double, more="gridSize(2)"); if (allocated(error)) return - call check(error, gridSize(3), 4.0_c_double, more="gridSize(3)"); if (allocated(error)) return + call check(error, gridSize(1), 1.0_ReKi, more="gridSize(1)"); if (allocated(error)) return + call check(error, gridSize(2), 2.0_ReKi, more="gridSize(2)"); if (allocated(error)) return + call check(error, gridSize(3), 4.0_ReKi, more="gridSize(3)"); if (allocated(error)) return ! Check lower bounds - call check(error, bounds(1,1), 6.5_c_double, more="origin(1)"); if (allocated(error)) return - call check(error, bounds(1,2), 6.5_c_double, more="origin(2)"); if (allocated(error)) return - call check(error, bounds(1,3), 6.5_c_double, more="origin(3)"); if (allocated(error)) return + call check(error, bounds(1,1), 6.5_ReKi, more="origin(1)"); if (allocated(error)) return + call check(error, bounds(1,2), 6.5_ReKi, more="origin(2)"); if (allocated(error)) return + call check(error, bounds(1,3), 6.5_ReKi, more="origin(3)"); if (allocated(error)) return ! Check upper bounds - call check(error, bounds(2,1), 7.5_c_double, more="ub(1)"); if (allocated(error)) return - call check(error, bounds(2,2), 8.5_c_double, more="ub(2)"); if (allocated(error)) return - call check(error, bounds(2,3), 10.5_c_double, more="ub(3)"); if (allocated(error)) return + call check(error, bounds(2,1), 7.5_ReKi, more="ub(1)"); if (allocated(error)) return + call check(error, bounds(2,2), 8.5_ReKi, more="ub(2)"); if (allocated(error)) return + call check(error, bounds(2,3), 10.5_ReKi, more="ub(3)"); if (allocated(error)) return ! Display grid properties print*, "dir = ", trim(DirPath) @@ -191,16 +211,20 @@ subroutine AMReX_test_amrex_find_subvols_1(error) integer(IntKi), parameter :: NumSteps = 5 character(*), parameter :: StartIndex = "00000" - integer(IntKi) :: FirstIndex, IndexDelta - integer(IntKi) :: ErrStat + integer(IntKi), allocatable :: DirIndices(:) + integer(IntKi), parameter :: Expected(0:NumSteps-1) = [0, 6, 12, 18, 24] + integer(IntKi) :: ErrStat, i character(ErrMsgLen) :: ErrMsg call amrex_find_subvols(DirPath, SubVol, DT, NumSteps, StartIndex, & - FirstIndex, IndexDelta, ErrStat, ErrMsg) + DirIndices, ErrStat, ErrMsg) call check(error, ErrStat, ErrID_None, more="amrex_find_subvols: "//trim(ErrMsg)); if (allocated(error)) return - call check(error, FirstIndex, 0); if (allocated(error)) return - call check(error, IndexDelta, 6); if (allocated(error)) return + call check(error, lbound(DirIndices,1), 0, more="lbound"); if (allocated(error)) return + call check(error, ubound(DirIndices,1), NumSteps-1, more="ubound"); if (allocated(error)) return + do i = 0, NumSteps-1 + call check(error, DirIndices(i), Expected(i), more="step "//trim(Num2LStr(i))); if (allocated(error)) return + end do end subroutine @@ -214,16 +238,18 @@ subroutine AMReX_test_amrex_find_subvols_2(error) integer(IntKi), parameter :: NumSteps = 3 character(*), parameter :: StartIndex = "00016" - integer(IntKi) :: FirstIndex, IndexDelta - integer(IntKi) :: ErrStat + integer(IntKi), allocatable :: DirIndices(:) + integer(IntKi), parameter :: Expected(0:NumSteps-1) = [16, 24, 32] + integer(IntKi) :: ErrStat, i character(ErrMsgLen) :: ErrMsg call amrex_find_subvols(DirPath, SubVol, DT, NumSteps, StartIndex, & - FirstIndex, IndexDelta, ErrStat, ErrMsg) + DirIndices, ErrStat, ErrMsg) call check(error, ErrStat, ErrID_None, more="amrex_find_subvols: "//trim(ErrMsg)); if (allocated(error)) return - call check(error, FirstIndex, 16); if (allocated(error)) return - call check(error, IndexDelta, 8); if (allocated(error)) return + do i = 0, NumSteps-1 + call check(error, DirIndices(i), Expected(i), more="step "//trim(Num2LStr(i))); if (allocated(error)) return + end do end subroutine @@ -238,16 +264,427 @@ subroutine AMReX_test_amrex_find_subvols_3(error) integer(IntKi), parameter :: NumSteps = 3 character(*), parameter :: StartIndex = "00006" - integer(IntKi) :: FirstIndex, IndexDelta + integer(IntKi), allocatable :: DirIndices(:) + integer(IntKi), parameter :: Expected(0:NumSteps-1) = [6, 18, 30] + integer(IntKi) :: ErrStat, i + character(ErrMsgLen) :: ErrMsg + + call amrex_find_subvols(DirPath, SubVol, DT, NumSteps, StartIndex, & + DirIndices, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="amrex_find_subvols: "//trim(ErrMsg)); if (allocated(error)) return + + do i = 0, NumSteps-1 + call check(error, DirIndices(i), Expected(i), more="step "//trim(Num2LStr(i))); if (allocated(error)) return + end do + + end subroutine + + ! Sub-volume directories written with a varying solver time step: the output interval is a + ! uniform 0.1 s but the step index stride changes from 4 to 2 part way through, as happens when + ! AMR-Wind transitions from time.initial_dt to fixed_dt. Header times carry the floating-point + ! drift observed in real data. Before directories were matched by time this failed with + ! "inconsistent delta between indices '00008' and '00010'". + subroutine AMReX_test_find_subvols_variable_stride(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data/subvolvardt" + integer(IntKi), parameter :: SubVol = 0 + real(DbKi), parameter :: DT = 0.1_DbKi + integer(IntKi), parameter :: NumSteps = 5 + character(*), parameter :: StartIndex = "00000" + + integer(IntKi), allocatable :: DirIndices(:) + integer(IntKi), parameter :: Expected(0:NumSteps-1) = [0, 4, 8, 10, 12] + integer(IntKi) :: ErrStat, i + character(ErrMsgLen) :: ErrMsg + + call amrex_find_subvols(DirPath, SubVol, DT, NumSteps, StartIndex, & + DirIndices, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="amrex_find_subvols: "//trim(ErrMsg)); if (allocated(error)) return + + call check(error, lbound(DirIndices,1), 0, more="lbound"); if (allocated(error)) return + do i = 0, NumSteps-1 + call check(error, DirIndices(i), Expected(i), more="step "//trim(Num2LStr(i))); if (allocated(error)) return + end do + + end subroutine + + ! Asking for one step more than exists must name the missing step, not report a bare count + subroutine AMReX_test_find_subvols_missing_step(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data/subvolvardt" + integer(IntKi), parameter :: SubVol = 0 + real(DbKi), parameter :: DT = 0.1_DbKi + integer(IntKi), parameter :: NumSteps = 6 + character(*), parameter :: StartIndex = "00000" + + integer(IntKi), allocatable :: DirIndices(:) + integer(IntKi) :: ErrStat + character(ErrMsgLen) :: ErrMsg + + call amrex_find_subvols(DirPath, SubVol, DT, NumSteps, StartIndex, & + DirIndices, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_Fatal, more="expected a fatal error"); if (allocated(error)) return + call check(error, index(ErrMsg, "time step 5") > 0, .true., & + more="message should name the missing step: "//trim(ErrMsg)); if (allocated(error)) return + + end subroutine + + ! Two directories whose header times land on the same step is ambiguous and must be rejected + subroutine AMReX_test_find_subvols_duplicate_step(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data/subvolvardt" + integer(IntKi), parameter :: SubVol = 1 + real(DbKi), parameter :: DT = 0.1_DbKi + integer(IntKi), parameter :: NumSteps = 2 + character(*), parameter :: StartIndex = "00000" + + integer(IntKi), allocatable :: DirIndices(:) + integer(IntKi) :: ErrStat + character(ErrMsgLen) :: ErrMsg + + call amrex_find_subvols(DirPath, SubVol, DT, NumSteps, StartIndex, & + DirIndices, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_Fatal, more="expected a fatal error"); if (allocated(error)) return + call check(error, index(ErrMsg, "claim time step 1") > 0, .true., & + more="message should report the duplicated step: "//trim(ErrMsg)); if (allocated(error)) return + + end subroutine + + ! A directory sitting off the step grid by more than the tolerance should be reported as a + ! near miss, so the user is pointed at the file rather than left with a bare count mismatch + subroutine AMReX_test_find_subvols_out_of_tolerance(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data/subvolvardt" + integer(IntKi), parameter :: SubVol = 2 + real(DbKi), parameter :: DT = 0.1_DbKi + integer(IntKi), parameter :: NumSteps = 2 + character(*), parameter :: StartIndex = "00000" + + integer(IntKi), allocatable :: DirIndices(:) + integer(IntKi) :: ErrStat + character(ErrMsgLen) :: ErrMsg + + call amrex_find_subvols(DirPath, SubVol, DT, NumSteps, StartIndex, & + DirIndices, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_Fatal, more="expected a fatal error"); if (allocated(error)) return + call check(error, index(ErrMsg, "outside") > 0, .true., & + more="message should flag the near miss: "//trim(ErrMsg)); if (allocated(error)) return + + end subroutine + + ! Directory indices that grow past the width of DirStartIndex. AMReX pads to a minimum width and + ! widens beyond it, so a six-digit index follows a five-digit one. Comparing the suffixes as text + ! silently drops those directories, because "100000" sorts before "99998". + subroutine AMReX_test_find_subvols_wide_index(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data/subvolwide" + integer(IntKi), parameter :: SubVol = 0 + real(DbKi), parameter :: DT = 0.1_DbKi + integer(IntKi), parameter :: NumSteps = 3 + character(*), parameter :: StartIndex = "99998" + + integer(IntKi), allocatable :: DirIndices(:) + integer(IntKi), parameter :: Expected(0:NumSteps-1) = [99998, 100000, 100002] + integer(IntKi) :: ErrStat, i + character(ErrMsgLen) :: ErrMsg + + call amrex_find_subvols(DirPath, SubVol, DT, NumSteps, StartIndex, & + DirIndices, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="amrex_find_subvols: "//trim(ErrMsg)); if (allocated(error)) return + + do i = 0, NumSteps-1 + call check(error, DirIndices(i), Expected(i), more="step "//trim(Num2LStr(i))); if (allocated(error)) return + end do + + end subroutine + + ! --------------------------------------------------------------------------------------- + ! ReadWindAMReX: resolving a time step to a directory + ! + ! These cover the contract between the index table built during initialization and the + ! read itself: the table is indexed by the 0-based FAST.Farm time step, and the directory + ! suffix is zero-padded to at least DirIndexLen characters and widened beyond it when the + ! index needs more digits. + ! --------------------------------------------------------------------------------------- + + ! Step n must resolve to the directory the table names for it, not to a computed index. + ! Table [0, 6, 12] means step 1 is directory 00006; read it and compare against reading + ! that directory directly. + subroutine AMReX_test_ReadWindAMReX_lookup(error) + type(error_type), allocatable, intent(out) :: error + type(AWAE_ParameterType) :: p + real(SiKi), allocatable :: viaTable(:,:,:,:), direct(:,:,:,:) + integer(IntKi) :: ErrStat + character(ErrMsgLen) :: ErrMsg + + ! Sub-volume 0 of the subvolmultiple fixture set is a 3x4x5 grid + allocate(viaTable(3,3,4,5)); allocate(direct(3,3,4,5)) + + p%WindFilePath = "data/subvolmultiple" + p%DirIndexLen = 5 + allocate(p%DirIndexLow(0:2)) + p%DirIndexLow = [0, 6, 12] + + call ReadWindAMReX(0, 1, p, viaTable, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="ReadWindAMReX: "//trim(ErrMsg)); if (allocated(error)) return + + call amrex_read_data("data/subvolmultiple_0_00006", direct, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="amrex_read_data: "//trim(ErrMsg)); if (allocated(error)) return + + call check(error, all(viaTable == direct), .true., & + more="step 1 did not resolve to directory 00006"); if (allocated(error)) return + + end subroutine + + ! A step outside the table is a hard error. It must never be clamped to the nearest + ! entry, which would silently freeze the inflow for the rest of the run. + subroutine AMReX_test_ReadWindAMReX_out_of_range(error) + type(error_type), allocatable, intent(out) :: error + type(AWAE_ParameterType) :: p + real(SiKi), allocatable :: dat(:,:,:,:) + integer(IntKi) :: ErrStat + character(ErrMsgLen) :: ErrMsg + + allocate(dat(3,3,4,5)) + p%WindFilePath = "data/subvolmultiple" + p%DirIndexLen = 5 + allocate(p%DirIndexLow(0:2)) + p%DirIndexLow = [0, 6, 12] + + call ReadWindAMReX(0, 3, p, dat, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_Fatal, more="expected a fatal error"); if (allocated(error)) return + call check(error, index(ErrMsg, "only steps 0 through 2") > 0, .true., & + more="message should report the available range: "//trim(ErrMsg)); if (allocated(error)) return + + end subroutine + + ! An index needing more digits than DirIndexLen must widen the field rather than overflow + ! it. The subvolwide fixtures run 99998 -> 100000 -> 100002 with a five-character start + ! index, so step 1 can only be read if the suffix widened to six digits; a fixed-width + ! write would have produced '*****' and failed to open anything. + subroutine AMReX_test_ReadWindAMReX_wide_index(error) + type(error_type), allocatable, intent(out) :: error + type(AWAE_ParameterType) :: p + real(SiKi), allocatable :: viaTable(:,:,:,:), direct(:,:,:,:) integer(IntKi) :: ErrStat character(ErrMsgLen) :: ErrMsg + allocate(viaTable(3,3,4,5)); allocate(direct(3,3,4,5)) + + p%WindFilePath = "data/subvolwide" + p%DirIndexLen = 5 + allocate(p%DirIndexLow(0:2)) + p%DirIndexLow = [99998, 100000, 100002] + + call ReadWindAMReX(0, 1, p, viaTable, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="ReadWindAMReX: "//trim(ErrMsg)); if (allocated(error)) return + + call amrex_read_data("data/subvolwide_0_100000", direct, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="amrex_read_data: "//trim(ErrMsg)); if (allocated(error)) return + + call check(error, all(viaTable == direct), .true., & + more="six-digit index did not resolve to directory 100000"); if (allocated(error)) return + + end subroutine + + ! If initialization never populated the table, say so instead of reading whatever a + ! zero-filled lookup would point at. + subroutine AMReX_test_ReadWindAMReX_no_table(error) + type(error_type), allocatable, intent(out) :: error + type(AWAE_ParameterType) :: p + real(SiKi), allocatable :: dat(:,:,:,:) + integer(IntKi) :: ErrStat + character(ErrMsgLen) :: ErrMsg + + allocate(dat(3,3,4,5)) + p%WindFilePath = "data/subvolmultiple" + p%DirIndexLen = 5 + + call ReadWindAMReX(0, 0, p, dat, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_Fatal, more="expected a fatal error"); if (allocated(error)) return + call check(error, index(ErrMsg, "never populated") > 0, .true., & + more="message should name the cause: "//trim(ErrMsg)); if (allocated(error)) return + + end subroutine + + ! --------------------------------------------------------------------------------------- + ! Scan robustness and the fast header path + ! --------------------------------------------------------------------------------------- + + ! A leftover directory from an earlier run can carry a time far past the window on a LOWER + ! index than valid data (different time step, same index base). The ascending-index walk + ! must not stop there while steps are still unclaimed: 00002 has t = 5.0 but 00004 and + ! 00008 hold steps 1 and 2. + subroutine AMReX_test_find_subvols_stale_beyond_window(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data/subvolstale" + integer(IntKi), parameter :: SubVol = 0 + real(DbKi), parameter :: DT = 0.1_DbKi + integer(IntKi), parameter :: NumSteps = 3 + character(*), parameter :: StartIndex = "00000" + + integer(IntKi), allocatable :: DirIndices(:) + integer(IntKi), parameter :: Expected(0:NumSteps-1) = [0, 4, 8] + integer(IntKi) :: ErrStat, i + character(ErrMsgLen) :: ErrMsg + + call amrex_find_subvols(DirPath, SubVol, DT, NumSteps, StartIndex, & + DirIndices, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="amrex_find_subvols: "//trim(ErrMsg)); if (allocated(error)) return + do i = 0, NumSteps-1 + call check(error, DirIndices(i), Expected(i), more="step "//trim(Num2LStr(i))); if (allocated(error)) return + end do + + end subroutine + + ! Real Kynema/AMR-Wind sub-volume output layout: one FAB tiling the domain box, three masked + ! velocity components, headers exactly as written by the sampler. Unlike the hand-built + ! subvolmultiple fixtures, the Header domain box equals the Cell_H box union here, so the + ! text-parse fast path engages. The FAB payloads were replaced by a known field + ! (u = 1+i, v = 2+j, w = 3+k in 0-based cell indices) so the read can be checked exactly. + subroutine AMReX_test_find_subvols_real_tiled(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data/subvoltiled" + integer(IntKi), parameter :: SubVol = 1 + real(DbKi), parameter :: DT = 0.1_DbKi + integer(IntKi), parameter :: NumSteps = 3 + character(*), parameter :: StartIndex = "31220" + + integer(IntKi), allocatable :: DirIndices(:) + integer(IntKi), parameter :: Expected(0:NumSteps-1) = [31220, 31224, 31228] + integer(IntKi) :: ErrStat, i + character(ErrMsgLen) :: ErrMsg + + call amrex_find_subvols(DirPath, SubVol, DT, NumSteps, StartIndex, & + DirIndices, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="amrex_find_subvols: "//trim(ErrMsg)); if (allocated(error)) return + do i = 0, NumSteps-1 + call check(error, DirIndices(i), Expected(i), more="step "//trim(Num2LStr(i))); if (allocated(error)) return + end do + + end subroutine + + subroutine AMReX_test_read_real_tiled(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data/subvoltiled_1_31220" + real(SiKi), allocatable :: data(:,:,:,:) + integer(IntKi) :: ErrStat, dims(3) + character(ErrMsgLen) :: ErrMsg + real(DbKi) :: time + real(ReKi) :: origin(3), gridSpacing(3) + + call amrex_read_header(DirPath, time, dims, gridSpacing, origin, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="amrex_read_header: "//trim(ErrMsg)); if (allocated(error)) return + call check(error, dims(1), 18, more="dims(1)"); if (allocated(error)) return + call check(error, dims(2), 18, more="dims(2)"); if (allocated(error)) return + call check(error, dims(3), 19, more="dims(3)"); if (allocated(error)) return + call check(error, gridSpacing(1), 8.0_ReKi, thr=1.0e-6_ReKi, more="dx"); if (allocated(error)) return + + allocate(data(3, dims(1), dims(2), dims(3))) + call amrex_read_data(DirPath, data, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="amrex_read_data: "//trim(ErrMsg)); if (allocated(error)) return + ! Known payload: u = 1+i, v = 2+j, w = 3+k with 0-based cell indices, checked at the corners + ! and an interior point. This also pins the component order and the x-fastest layout. + call check(error, data(1,1,1,1), 1.0_SiKi, thr=1.0e-6_SiKi, more="u(0,0,0)"); if (allocated(error)) return + call check(error, data(2,1,1,1), 2.0_SiKi, thr=1.0e-6_SiKi, more="v(0,0,0)"); if (allocated(error)) return + call check(error, data(3,1,1,1), 3.0_SiKi, thr=1.0e-6_SiKi, more="w(0,0,0)"); if (allocated(error)) return + call check(error, data(1,18,1,1), 18.0_SiKi, thr=1.0e-6_SiKi, more="u(17,0,0)"); if (allocated(error)) return + call check(error, data(2,1,18,1), 19.0_SiKi, thr=1.0e-6_SiKi, more="v(0,17,0)"); if (allocated(error)) return + call check(error, data(3,1,1,19), 21.0_SiKi, thr=1.0e-6_SiKi, more="w(0,0,18)"); if (allocated(error)) return + call check(error, data(1,5,7,9), 5.0_SiKi, thr=1.0e-6_SiKi, more="u(4,6,8)"); if (allocated(error)) return + call check(error, data(2,5,7,9), 8.0_SiKi, thr=1.0e-6_SiKi, more="v(4,6,8)"); if (allocated(error)) return + call check(error, data(3,5,7,9), 11.0_SiKi, thr=1.0e-6_SiKi, more="w(4,6,8)"); if (allocated(error)) return + + end subroutine + + ! The fast path is only trusted after the text parse agrees with amrex_read_header on the + ! starting directory. On real sub-volume output it must agree exactly. + subroutine AMReX_test_header_text_agrees_on_tiled(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data/subvoltiled_1_31220" + logical :: ok + integer(IntKi) :: ErrStat, dimsA(3), dimsT(3), i + character(ErrMsgLen) :: ErrMsg + real(DbKi) :: timeA, timeT + real(ReKi) :: originA(3), dxA(3), originT(3), dxT(3) + + call amrex_read_header(DirPath, timeA, dimsA, dxA, originA, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="amrex_read_header: "//trim(ErrMsg)); if (allocated(error)) return + call amrex_parse_header_text(DirPath, ok, timeT, dimsT, dxT, originT) + call check(error, ok, .true., more="text parse should succeed on real output"); if (allocated(error)) return + call check(error, timeT, timeA, thr=1.0e-9_DbKi, more="time"); if (allocated(error)) return + do i = 1, 3 + call check(error, dimsT(i), dimsA(i), more="dims"); if (allocated(error)) return + call check(error, dxT(i), dxA(i), thr=1.0e-6_ReKi, more="dx"); if (allocated(error)) return + call check(error, originT(i), originA(i), thr=1.0e-3_ReKi, more="origin"); if (allocated(error)) return + end do + + end subroutine + + ! The hand-built fixtures store a small box array inside a much larger domain box, so the + ! text parse (domain) and amrex_read_header (box union) disagree on dims. That disagreement + ! is exactly what makes the search fall back to the slow path for them -- record it. + subroutine AMReX_test_header_text_disagrees_on_subset(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data/subvolmultiple_0_00006" + logical :: ok + integer(IntKi) :: ErrStat, dimsA(3), dimsT(3) + character(ErrMsgLen) :: ErrMsg + real(DbKi) :: timeA, timeT + real(ReKi) :: originA(3), dxA(3), originT(3), dxT(3) + + call amrex_read_header(DirPath, timeA, dimsA, dxA, originA, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_None, more="amrex_read_header: "//trim(ErrMsg)); if (allocated(error)) return + call amrex_parse_header_text(DirPath, ok, timeT, dimsT, dxT, originT) + call check(error, ok, .true., more="text parse should succeed"); if (allocated(error)) return + call check(error, timeT, timeA, thr=1.0e-9_DbKi, more="time agrees"); if (allocated(error)) return + call check(error, dimsT(1), 256, more="text dims are the domain box"); if (allocated(error)) return + call check(error, dimsA(1), 3, more="reader dims are the box union"); if (allocated(error)) return + + end subroutine + + ! The scan lists the parent directory and matches its entries against the prefix. The entry + ! paths the iterator hands back are not textually the prefix the caller passed -- a prefix with + ! no directory of its own comes back as "./name", and a redundant separator is normalized away + ! -- so the comparison has to be made on the final path component. Matching the whole path + ! instead finds nothing beyond the starting directory even though every directory is present. + subroutine AMReX_test_find_subvols_prefix_form(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data//subvolmultiple" + integer(IntKi), parameter :: SubVol = 0 + real(DbKi), parameter :: DT = 0.6_DbKi + integer(IntKi), parameter :: NumSteps = 5 + character(*), parameter :: StartIndex = "00000" + + integer(IntKi), allocatable :: DirIndices(:) + integer(IntKi), parameter :: Expected(0:NumSteps-1) = [0, 6, 12, 18, 24] + integer(IntKi) :: ErrStat, i + character(ErrMsgLen) :: ErrMsg + call amrex_find_subvols(DirPath, SubVol, DT, NumSteps, StartIndex, & - FirstIndex, IndexDelta, ErrStat, ErrMsg) + DirIndices, ErrStat, ErrMsg) call check(error, ErrStat, ErrID_None, more="amrex_find_subvols: "//trim(ErrMsg)); if (allocated(error)) return + do i = 0, NumSteps-1 + call check(error, DirIndices(i), Expected(i), more="step "//trim(Num2LStr(i))); if (allocated(error)) return + end do + + end subroutine + + ! amrex_read_data writes the destination by grid index, so a plotfile describing a different + ! grid must be rejected rather than written past the end of the caller's array. The fixture + ! grid is 3x4x5; ask for it into a 3x4x4 array. + subroutine AMReX_test_read_data_wrong_dims(error) + type(error_type), allocatable, intent(out) :: error + character(*), parameter :: DirPath = "data/subvolmultiple_0_00006" + real(SiKi), allocatable :: dat(:,:,:,:) + integer(IntKi) :: ErrStat + character(ErrMsgLen) :: ErrMsg - call check(error, FirstIndex, 6, more="FirstIndex"); if (allocated(error)) return - call check(error, IndexDelta, 12, more="IndexDelta"); if (allocated(error)) return + allocate(dat(3,3,4,4)) + call amrex_read_data(DirPath, dat, ErrStat, ErrMsg) + call check(error, ErrStat, ErrID_Fatal, more="expected a fatal error"); if (allocated(error)) return + call check(error, index(ErrMsg, "do not match") > 0, .true., & + more="message should report the dimension mismatch: "//trim(ErrMsg)); if (allocated(error)) return end subroutine diff --git a/modules/wakedynamics/src/WakeDynamics.f90 b/modules/wakedynamics/src/WakeDynamics.f90 index 05b357088e..613b4999e1 100644 --- a/modules/wakedynamics/src/WakeDynamics.f90 +++ b/modules/wakedynamics/src/WakeDynamics.f90 @@ -45,8 +45,8 @@ module WakeDynamics public :: WD_WritePlaneOutputs ! Routine for IO Operation public :: WD_CalcConstrStateResidual ! Tight coupling routine for returning the constraint state residual - public :: WD_TEST_Axi2Cart - public :: WD_TEST_AddVelocityCurl + public :: AddVelocityCurl ! Exposed for unit testing + public :: Axisymmetric2CartesianVel ! Exposed for unit testing contains function WD_Interp ( yVal, xArr, yArr ) @@ -426,6 +426,7 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut !............................................................................................ p%TurbNum = InitInp%TurbNum p%DT_low = interval + p%LowResBounds = InitInp%LowResBounds ! Parameters from input file p%Mod_Wake = InitInp%InputFileData%Mod_Wake p%MaxNumPlanes = InitInp%MaxNumPlanes @@ -500,7 +501,9 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut allocate( u%Ct_azavg ( 0:p%NumRadii-1 ),stat=errStat2); if (Failed0('u%Ct_azavg.')) return; allocate( u%Cq_azavg ( 0:p%NumRadii-1 ),stat=errStat2); if (Failed0('u%Cq_azavg.')) return; if (errStat /= ErrID_None) return - + u%V_plane = 0.0_ReKi + u%Ct_azavg = 0.0_ReKi + u%Cq_azavg = 0.0_ReKi @@ -548,6 +551,8 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut xd%Vx_wake = 0.0_ReKi xd%Vr_wake = 0.0_ReKi xd%Vx_wake2 = 0.0_ReKi + xd%Vy_wake2 = 0.0_ReKi + xd%Vz_wake2 = 0.0_ReKi xd%V_plane_filt = 0.0_ReKi xd%Vx_wind_disk_filt = 0.0_ReKi xd%TI_amb_filt = 0.0_ReKi @@ -556,6 +561,7 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut xd%Ct_azavg_filt = 0.0_ReKi xd%Cq_azavg_filt = 0.0_ReKi OtherState%firstPass = .true. + OtherState%MaxPlanesWarned = .false. ! miscvars to avoid the allocation per timestep ! Cartesian eddy viscosity (allocated even for polar if plane outputs are requested) @@ -574,12 +580,15 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut allocate ( m%vt_tot (0:p%NumRadii-1,0:p%MaxNumPlanes-1 ) , STAT=ErrStat2 ); if (Failed0('m%vt_tot.')) return; allocate ( m%vt_amb (0:p%NumRadii-1,0:p%MaxNumPlanes-1 ) , STAT=ErrStat2 ); if (Failed0('m%vt_amb.')) return; allocate ( m%vt_shr (0:p%NumRadii-1,0:p%MaxNumPlanes-1 ) , STAT=ErrStat2 ); if (Failed0('m%vt_shr.')) return; + m%dvtdr = 0.0_ReKi + m%vt_tot = 0.0_ReKi + m%vt_amb = 0.0_ReKi + m%vt_shr = 0.0_ReKi else if (p%Mod_Wake == Mod_Wake_Cartesian .or. p%Mod_Wake == Mod_Wake_Curl) then allocate ( m%nu_dvx_dy(-p%NumRadii+1:p%NumRadii-1,-p%NumRadii+1:p%NumRadii-1), STAT=ErrStat2 ); if (Failed0('m%nu_dvx_dy.')) return; allocate ( m%nu_dvx_dz(-p%NumRadii+1:p%NumRadii-1,-p%NumRadii+1:p%NumRadii-1), STAT=ErrStat2 ); if (Failed0('m%nu_dvx_dz.')) return; allocate ( m%dnuvx_dy (-p%NumRadii+1:p%NumRadii-1,-p%NumRadii+1:p%NumRadii-1), STAT=ErrStat2 ); if (Failed0('m%dnuvx_dy.' )) return; allocate ( m%dnuvx_dz (-p%NumRadii+1:p%NumRadii-1,-p%NumRadii+1:p%NumRadii-1), STAT=ErrStat2 ); if (Failed0('m%dnuvx_dz.' )) return; - if (errStat /= ErrID_None) return m%nu_dvx_dy = 0.0_ReKi m%nu_dvx_dz = 0.0_ReKi m%dnuvx_dy = 0.0_ReKi @@ -595,10 +604,18 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut allocate ( m%d(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%d.')) return; allocate ( m%r_wake(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%r_wake.' )) return; allocate ( m%Vx_high(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%Vx_high.' )) return; - allocate ( m%Vt_wake(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%Vx_high.' )) return; + allocate ( m%Vt_wake(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%Vt_wake.' )) return; allocate ( m%Vx_polar(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%Vx_polar.')) return; + + m%a = 0.0_ReKi + m%b = 0.0_ReKi + m%c = 0.0_ReKi + m%d = 0.0_ReKi + m%r_wake = 0.0_ReKi + m%Vx_high = 0.0_ReKi m%Vx_polar = 0.0_ReKi - m%Vt_wake = 0.0_ReKi + m%Vt_wake = 0.0_ReKi + !............................................................................................ ! Define initialization output here !............................................................................................ @@ -715,6 +732,9 @@ subroutine WD_UpdateStates( t, n, u, p, x, xd, z, OtherState, m, errStat, errMsg integer(intKi) :: i,j, maxPln integer(intKi) :: iy, iz ! indices on y and z real(ReKi) :: vt_min ! Minimum Eddy viscosity + integer(IntKi) :: oobIdx(0:p%MaxNumPlanes) ! One extra slot: allows all p%MaxNumPlanes planes can be simultaneously out-of-bounds (corner case) + integer(IntKi) :: nOOB, iOOB, jOOB + logical :: merged errStat = ErrID_None errMsg = "" @@ -950,90 +970,153 @@ subroutine WD_UpdateStates( t, n, u, p, x, xd, z, OtherState, m, errStat, errMsg endif endif - !Used for debugging: write(51,'(I5,100(1x,ES10.2E2))') n, xd%x_plane(n), xd%x_plane(n)/xd%D_rotor_filt(n), xd%Vx_wind_disk_filt(n) + xd%Vx_wake(:,n), xd%Vr_wake(:,n) + + ! -------------------------------------------------------------------------------- + ! Drop planes that exit the allocated memory buffer + ! -------------------------------------------------------------------------------- xd%NumPlanes = xd%NumPlanes + 1.0 if ( NINT(xd%NumPlanes) > p%MaxNumPlanes ) then xd%NumPlanes = real(p%MaxNumPlanes,ReKi) - call SetErrStat(ErrID_Warn, ' The number of wake planes of turbine '//trim(num2lstr(p%TurbNum))//' exceeded the allowed number ('//trim(num2lstr(p%MaxNumPlanes))//'). Excess plane(s) removed. ', errStat, errMsg, RoutineName) - if (errStat >= AbortErrLev) then - call Cleanup() - return + if (.not. OtherState%MaxPlanesWarned) then + call SetErrStat(ErrID_Warn, ' The number of wake planes of turbine '//trim(num2lstr(p%TurbNum))//' exceeded the allowed number ('//trim(num2lstr(p%MaxNumPlanes))//'). Excess plane(s) removed. ', errStat, errMsg, RoutineName) + OtherState%MaxPlanesWarned = .true. end if end if if ( NINT(xd%NumPlanes) < 2 ) then ! Check just in case following implementation plan; however, this should never happen. Consider removing in the future. call SetErrStat(ErrID_Fatal, ' The number of wake planes of turbine '//trim(num2lstr(p%TurbNum))//' has dropped below 2. ', errStat, errMsg, RoutineName) + call Cleanup() + return end if + ! -------------------------------------------------------------------------------- + ! Merge any out-of-bounds planes within 2*dr of each other (non-sequential) + ! -------------------------------------------------------------------------------- maxPln = NINT(xd%NumPlanes) - 1 - do i=maxPln,0,-1 - - ! if a plane is beyond the buffer, simply drop it and all following planes (it should only be the last plane that gets dropped) - if ( xd%x_plane(i) > p%x_Buff ) then - xd%NumPlanes = max( xd%NumPlanes - 1.0, 2.0 ) ! Plane indexing includes 0, hence the -1.0 - cycle - endif - - ! If a plane overtakes another plane, merge the planes by averaging, then shift all remaining planes forward. - if ( i+1 < NINT(xd%NumPlanes)) then ! don't overstep bounds with i+1 indexing - if (xd%x_plane(i) >= xd%x_plane(i+1) ) then - - call SetErrStat(ErrID_Warn, ' Turbine '//trim(num2lstr(p%TurbNum))//' wake plane '//trim(num2lstr(i))// & - ' (x_plane='//trim(num2lstr(xd%x_plane(i)))//') has overtaken wake plane '//trim(num2lstr(i+1))// & - ' (x_plane='//trim(num2lstr(xd%x_plane(i+1)))// & - '). Merging planes by averaging. Reduce f_c to prevent planes from passing each other. ', errStat, errMsg, RoutineName) - if (errStat >= AbortErrLev) then - call Cleanup() - return - end if - - ! Merge the i and i+1 plane by averaging them together - xd%Vx_wind_disk_filt(i) = (xd%Vx_wind_disk_filt(i) + xd%Vx_wind_disk_filt(i+1)) / 2.0_ReKi - xd%x_plane ( i) = (xd%x_plane ( i) + xd%x_plane ( i+1)) / 2.0_ReKi - xd%TI_amb_filt ( i) = (xd%TI_amb_filt ( i) + xd%TI_amb_filt ( i+1)) / 2.0_ReKi - xd%D_rotor_filt ( i) = (xd%D_rotor_filt ( i) + xd%D_rotor_filt ( i+1)) / 2.0_ReKi - xd%YawErr_filt ( i) = (xd%YawErr_filt ( i) + xd%YawErr_filt ( i+1)) / 2.0_ReKi - xd%p_plane ( :,i) = (xd%p_plane ( :,i) + xd%p_plane ( :,i+1)) / 2.0_ReKi - xd%xhat_plane ( :,i) = (xd%xhat_plane ( :,i) + xd%xhat_plane ( :,i+1)) / 2.0_ReKi - xd%xhat_plane ( :,i) = xd%xhat_plane ( :,i) / TwoNorm(xd%xhat_plane( :,i)) ! renormalize - xd%V_plane_filt ( :,i) = (xd%V_plane_filt ( :,i) + xd%V_plane_filt ( :,i+1)) / 2.0_ReKi - xd%Vx_wake ( :,i) = (xd%Vx_wake ( :,i) + xd%Vx_wake ( :,i+1)) / 2.0_ReKi - xd%Vr_wake ( :,i) = (xd%Vr_wake ( :,i) + xd%Vr_wake ( :,i+1)) / 2.0_ReKi - xd%Vx_wake2 (:,:,i) = (xd%Vx_wake2 (:,:,i) + xd%Vx_wake2 (:,:,i+1)) / 2.0_ReKi - xd%Vy_wake2 (:,:,i) = (xd%Vy_wake2 (:,:,i) + xd%Vy_wake2 (:,:,i+1)) / 2.0_ReKi - xd%Vz_wake2 (:,:,i) = (xd%Vz_wake2 (:,:,i) + xd%Vz_wake2 (:,:,i+1)) / 2.0_ReKi - - ! Since i and i+1 planes are now merged effectively dropping a plane, shift all planes that follow forward - do j = i+1,NINT(xd%NumPlanes)-2 ! NumPlanes includes 0 index plane, so last valid index is NumPlanes-1. - xd%Vx_wind_disk_filt(j) = xd%Vx_wind_disk_filt(j+1) - xd%x_plane ( j) = xd%x_plane ( j+1) - xd%TI_amb_filt ( j) = xd%TI_amb_filt ( j+1) - xd%D_rotor_filt ( j) = xd%D_rotor_filt ( j+1) - xd%YawErr_filt ( j) = xd%YawErr_filt ( j+1) - xd%p_plane ( :,j) = xd%p_plane ( :,j+1) - xd%xhat_plane ( :,j) = xd%xhat_plane ( :,j+1) - xd%V_plane_filt ( :,j) = xd%V_plane_filt ( :,j+1) - xd%Vx_wake ( :,j) = xd%Vx_wake ( :,j+1) - xd%Vr_wake ( :,j) = xd%Vr_wake ( :,j+1) - xd%Vx_wake2 (:,:,j) = xd%Vx_wake2 (:,:,j+1) - xd%Vy_wake2 (:,:,j) = xd%Vy_wake2 (:,:,j+1) - xd%Vz_wake2 (:,:,j) = xd%Vz_wake2 (:,:,j+1) - end do - - ! Now that we shifted the planes up, remove the last one - xd%NumPlanes = xd%NumPlanes - 1.0 + ! Collect indices of all OOB planes + nOOB = 0 + do i = 0, maxPln + if (PlaneOutOfBounds(xd%p_plane(:,i))) then + nOOB = nOOB + 1 + oobIdx(nOOB) = i + end if + end do + ! Check all OOB pairs for proximity; work backwards so shifts don't invalidate lower indices + iOOB = nOOB + do while (iOOB >= 2) + merged = .false. + do jOOB = iOOB - 1, 1, -1 + if (TwoNorm(xd%p_plane(:,oobIdx(iOOB)) - xd%p_plane(:,oobIdx(jOOB))) <= 2.0_ReKi * p%dr) then + call MergeWakePlanes(oobIdx(jOOB), oobIdx(iOOB)) + ! Remove entry iOOB and adjust indices above the dropped plane + call AdjustOobIndices(oobIdx, nOOB, iOOB) + ! nOOB has shrunk; clamp iOOB so it still refers to a valid, populated entry + iOOB = min(iOOB, nOOB) + merged = .true. + exit end if - end if + end do + if (.not. merged) iOOB = iOOB - 1 + end do + ! -------------------------------------------------------------------------------- + ! Drop individual planes that have gone beyond the buffer. + ! Each out-of-buffer plane is removed and higher-indexed planes are shifted down, + ! so this works even when the furthest-travelled plane is not the last index + ! (e.g. for drones). + ! -------------------------------------------------------------------------------- + i = NINT(xd%NumPlanes) - 1 + do while (i >= 0 .and. NINT(xd%NumPlanes) > 2) + if ( xd%x_plane(i) > p%x_Buff ) then + call ShiftWakePlanesDown(i) + xd%NumPlanes = max( xd%NumPlanes - 1.0, 2.0 ) + ! Don't decrement i: the plane that shifted into position i needs checking too + else + i = i - 1 + endif end do call Cleanup() contains + !> Merge two adjacent wake planes by averaging their states into iKeep, + !! then shift all planes above iDrop down by one and decrement NumPlanes. + !! iKeep is the plane that survives (receives the average), iDrop is removed. + !! Typically iKeep = min(iA,iB) and iDrop = max(iA,iB). + subroutine MergeWakePlanes(iKeep, iDrop) + integer(IntKi), intent(in) :: iKeep !< Index of plane to keep (receives averaged values) + integer(IntKi), intent(in) :: iDrop !< Index of plane to remove (shifted out) + integer(IntKi) :: j + + ! Average the two planes into iKeep + xd%Vx_wind_disk_filt(iKeep) = (xd%Vx_wind_disk_filt(iKeep) + xd%Vx_wind_disk_filt(iDrop)) / 2.0_ReKi + xd%x_plane ( iKeep) = (xd%x_plane ( iKeep) + xd%x_plane ( iDrop)) / 2.0_ReKi + xd%TI_amb_filt ( iKeep) = (xd%TI_amb_filt ( iKeep) + xd%TI_amb_filt ( iDrop)) / 2.0_ReKi + xd%D_rotor_filt ( iKeep) = (xd%D_rotor_filt ( iKeep) + xd%D_rotor_filt ( iDrop)) / 2.0_ReKi + xd%YawErr_filt ( iKeep) = (xd%YawErr_filt ( iKeep) + xd%YawErr_filt ( iDrop)) / 2.0_ReKi + xd%p_plane ( :,iKeep) = (xd%p_plane ( :,iKeep) + xd%p_plane ( :,iDrop)) / 2.0_ReKi + xd%xhat_plane ( :,iKeep) = (xd%xhat_plane ( :,iKeep) + xd%xhat_plane ( :,iDrop)) / 2.0_ReKi + xd%xhat_plane ( :,iKeep) = xd%xhat_plane ( :,iKeep) / TwoNorm(xd%xhat_plane(:,iKeep)) + xd%V_plane_filt ( :,iKeep) = (xd%V_plane_filt ( :,iKeep) + xd%V_plane_filt ( :,iDrop)) / 2.0_ReKi + xd%Vx_wake ( :,iKeep) = (xd%Vx_wake ( :,iKeep) + xd%Vx_wake ( :,iDrop)) / 2.0_ReKi + xd%Vr_wake ( :,iKeep) = (xd%Vr_wake ( :,iKeep) + xd%Vr_wake ( :,iDrop)) / 2.0_ReKi + xd%Vx_wake2 (:,:,iKeep) = (xd%Vx_wake2 (:,:,iKeep) + xd%Vx_wake2 (:,:,iDrop)) / 2.0_ReKi + xd%Vy_wake2 (:,:,iKeep) = (xd%Vy_wake2 (:,:,iKeep) + xd%Vy_wake2 (:,:,iDrop)) / 2.0_ReKi + xd%Vz_wake2 (:,:,iKeep) = (xd%Vz_wake2 (:,:,iKeep) + xd%Vz_wake2 (:,:,iDrop)) / 2.0_ReKi + + ! Shift all planes above iDrop down by one and decrement plane count + call ShiftWakePlanesDown(iDrop) + xd%NumPlanes = xd%NumPlanes - 1.0 + end subroutine MergeWakePlanes + + !> Shift all wake-plane state arrays above index iDrop down by one. + !! This removes the plane at iDrop; the caller must also decrement xd%NumPlanes. + subroutine ShiftWakePlanesDown(iDrop) + integer(IntKi), intent(in) :: iDrop !< Index of plane to remove + integer(IntKi) :: j + do j = iDrop, NINT(xd%NumPlanes)-2 + xd%Vx_wind_disk_filt(j) = xd%Vx_wind_disk_filt(j+1) + xd%x_plane ( j) = xd%x_plane ( j+1) + xd%TI_amb_filt ( j) = xd%TI_amb_filt ( j+1) + xd%D_rotor_filt ( j) = xd%D_rotor_filt ( j+1) + xd%YawErr_filt ( j) = xd%YawErr_filt ( j+1) + xd%p_plane ( :,j) = xd%p_plane ( :,j+1) + xd%xhat_plane ( :,j) = xd%xhat_plane ( :,j+1) + xd%V_plane_filt ( :,j) = xd%V_plane_filt ( :,j+1) + xd%Vx_wake ( :,j) = xd%Vx_wake ( :,j+1) + xd%Vr_wake ( :,j) = xd%Vr_wake ( :,j+1) + xd%Vx_wake2 (:,:,j) = xd%Vx_wake2 (:,:,j+1) + xd%Vy_wake2 (:,:,j) = xd%Vy_wake2 (:,:,j+1) + xd%Vz_wake2 (:,:,j) = xd%Vz_wake2 (:,:,j+1) + end do + end subroutine ShiftWakePlanesDown + + !> Check whether a wake plane center is outside the low-resolution domain bounds. + pure function PlaneOutOfBounds(p_pos) result(outOfBounds) + real(ReKi), intent(in) :: p_pos(3) !< Plane center position (XYZ) + logical :: outOfBounds + outOfBounds = any(p_pos < p%LowResBounds(:,1)) .or. any(p_pos > p%LowResBounds(:,2)) + end function PlaneOutOfBounds + + !> Remove entry iRemoved from oobIdx and decrement stored indices above the dropped plane. + subroutine AdjustOobIndices(oobIdx, nOOB, iRemoved) + integer(IntKi), intent(inout) :: oobIdx(0:), nOOB + integer(IntKi), intent(in) :: iRemoved + integer(IntKi) :: droppedPlane, k + droppedPlane = oobIdx(iRemoved) + do k = iRemoved, nOOB - 1 + oobIdx(k) = oobIdx(k+1) + end do + nOOB = nOOB - 1 + do k = 1, nOOB + if (oobIdx(k) > droppedPlane) oobIdx(k) = oobIdx(k) - 1 + end do + end subroutine AdjustOobIndices + subroutine updateVelocityPolar() integer(intKi) :: i,j real(ReKi) :: dx, absdx @@ -1123,8 +1206,8 @@ subroutine updateVelocityCartesian() call gradient_z(m%nu_dvx_dz, p%dr, m%dnuvx_dz ) ! Loop through all the points on the plane (y, z) - do iz = -p%NumRadii+2, p%NumRadii-2 - do iy = -p%NumRadii+2, p%NumRadii-2 + do iz = -p%NumRadii+1, p%NumRadii-1 + do iy = -p%NumRadii+1, p%NumRadii-1 ! Eddy viscosity term divTau = m%dnuvx_dy(iy,iz) + m%dnuvx_dz(iy,iz) @@ -1308,28 +1391,6 @@ subroutine AddSwirl(r, Vt_wake, y, z, Vy_curl, Vz_curl) end subroutine AddSwirl -!> Test the curled wake velocity curl function -subroutine WD_TEST_AddVelocityCurl() - - real(ReKi) :: Vy_curl(2,2)=0.0_ReKi - real(ReKi) :: Vz_curl(2,2)=0.0_ReKi - real(ReKi) :: y(2)=(/ 0., 2./) - real(ReKi) :: z(2)=(/-1.,1./) - real(ReKi) :: Gamma0 - - call AddVelocityCurl(Vx=10., yaw_angle=0.1, nVortex=100, R=63., psi_skew=0.2, & - y=y, z=z, Ct_avg=0.7, sigma_d=0.2, Vy_curl=Vy_curl, Vz_curl=Vz_curl, Gamma0=Gamma0) - - if (abs(Vy_curl(1,1)+0.217109)>1e-4) then - print*,'Test fail for vy' - !STOP - endif - if (abs(Vz_curl(2,2)+4.459746e-2)>1e-4) then - print*,'>>> Test fail for vz' - !STOP - endif -end subroutine - !> Weighted average of two angles @@ -1505,40 +1566,6 @@ function exp_safe(x) endif end function exp_safe -subroutine WD_TEST_Axi2Cart() - real(ReKi) :: r(4)=(/0.,1.,2.,3./) -! real(ReKi) :: y(4)=(/-1.,0.,1.5,2./) -! real(ReKi) :: z(5)=(/-2.5,-1.5,0.,1.5,2./) - real(ReKi) :: y(4)=(/0.,1. ,1.5, 2./) - real(ReKi) :: z(5)=(/0.,0.5,1. ,1.5,2./) - real(ReKi) :: Vr_axi(4) - real(ReKi) :: Vx_axi(4) - real(ReKi) :: Vx(4,5)=0.0_ReKi - real(ReKi) :: Vy(4,5)=0.0_ReKi - real(ReKi) :: Vz(4,5)=0.0_ReKi - integer :: i,j - real(ReKi) :: Vr, r_tmp - Vr_axi=4._ReKi*r - Vx_axi=3._ReKi*r - call Axisymmetric2CartesianVel(Vx_axi, Vr_axi, r, y, z, Vx, Vy, Vz) - - do i = 1,size(y) - do j = 1,size(z) - r_tmp = sqrt(y(i)**2+z(j)**2) - Vr = sqrt(Vy(i,j)**2 + Vz(i,j)**2) - if (abs(Vr-4*r_tmp)>1e-3) then - print*,'>>Error Axi2Cart Vr',Vr,4*r_tmp - STOP - endif - if (abs(Vx(i,j)-3*r_tmp)>1e-3) then - print*,'>>Error Axi2Cart Vx',Vx(i,j),3*r_tmp - STOP - endif - enddo - enddo -end subroutine - - !---------------------------------------------------------------------------------------------------------------------------------- !> Routine for computing outputs, used in both loose and tight coupling. @@ -1616,7 +1643,17 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) if ( p%OutAllPlanes ) then call mkdir(p%OutFileVTKDir) endif - + + ! set initial outputs for Cartesian/Curl wake models from the polar info. + if (p%Mod_Wake == Mod_Wake_Cartesian .or. p%Mod_Wake == Mod_Wake_Curl) then + call Axisymmetric2CartesianVx(y%Vx_wake(:,0), p%r, p%y, p%z, y%Vx_wake2(:,:,0)) + y%Vy_wake2(:,:,0) = 0.0_ReKi + y%Vz_wake2(:,:,0) = 0.0_ReKi + y%Vx_wake2(:,:,1) = y%Vx_wake2(:,:,0) + y%Vy_wake2(:,:,1) = 0.0_ReKi + y%Vz_wake2(:,:,1) = 0.0_ReKi + end if + else y%x_plane = xd%x_plane y%p_plane = xd%p_plane @@ -1631,13 +1668,14 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) y%NumPlanes = xd%NumPlanes end if - ! --- Linearly decay wake deficits in the buffer region based on distance + ! --- Linearly decay wake deficits in the buffer region based on distance. + ! Polar arrays are tapered first; for the Polar wake model these tapered values + ! are then converted to the Cartesian output below, so the taper carries through. do i = 0,maxPln - if ( xd%x_plane(i) > p%x_full ) then - ! Note: Clamp to zero just in case, but all wake planes that propagated past x_buff should have been removed. - ScBuff = max( ( p%x_buff - xd%x_plane(i) ) / p%d_buff , 0.0 ) - y%Vx_wake(:,i) = y%Vx_wake(:,i) * ScBuff - y%Vr_wake(:,i) = y%Vr_wake(:,i) * ScBuff + ScBuff = BufferScale(xd%x_plane(i)) + if ( ScBuff < 1.0_ReKi ) then + y%Vx_wake(:,i) = y%Vx_wake(:,i) * ScBuff + y%Vr_wake(:,i) = y%Vr_wake(:,i) * ScBuff end if end do @@ -1662,11 +1700,29 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) enddo endif else if (p%Mod_Wake == Mod_Wake_Cartesian .or. p%Mod_Wake == Mod_Wake_Curl) then - do i = 0, maxPln - y%Vx_wake2(:,:,i) = xd%Vx_wake2(:,:,i) - y%Vy_wake2(:,:,i) = xd%Vy_wake2(:,:,i) - y%Vz_wake2(:,:,i) = xd%Vz_wake2(:,:,i) - enddo + ! For Cartesian/Curl, the Cartesian output arrays are copied from the (untapered) + ! discrete state. Apply the buffer-region taper here so that the wake deficit + ! handed to AWAE fades linearly from x_Full to x_Buff. Do NOT modify xd%*_wake2: + ! tapering the state would compound the scaling at every time step. + ! Skip on firstPass: y%Vx_wake2 was already set from NearWakeCorrection above, + ! and xd%Vx_wake2 has not yet been initialized (still zero). + if (.not. OtherState%firstPass) then + do i = 0, maxPln + ScBuff = BufferScale(xd%x_plane(i)) + y%Vx_wake2(:,:,i) = xd%Vx_wake2(:,:,i) * ScBuff + y%Vy_wake2(:,:,i) = xd%Vy_wake2(:,:,i) * ScBuff + y%Vz_wake2(:,:,i) = xd%Vz_wake2(:,:,i) * ScBuff + enddo + ! Recompute Cartesian gradients from the tapered output field so that the WAT + ! gradient term in Calc_k_WAT and the dvx_dy/dz VTK diagnostics fade consistently + ! with the velocity deficit in the buffer region. + if ( p%WAT .or. p%OutAllPlanes ) then + do i = 0, maxPln + call gradient_y(y%Vx_wake2(:,:,i), p%dr, m%dvx_dy(:,:,i)) + call gradient_z(y%Vx_wake2(:,:,i), p%dr, m%dvx_dz(:,:,i)) + end do + endif + end if endif ! Curl or Polar ! --- WAT - Compute k_mt and add turbulence @@ -1675,6 +1731,22 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) end if contains + !> Linear buffer-region taper applied to the output wake fields. + !! Returns 1 for x_plane <= x_Full, fades linearly to 0 at x_Buff, and is + !! clamped to 0 beyond x_Buff (or whenever d_Buff <= 0, which avoids a + !! divide-by-zero when the user sets NumDBuff = 0). + pure function BufferScale(x_plane) result(ScBuff_loc) + real(ReKi), intent(in) :: x_plane + real(ReKi) :: ScBuff_loc + if ( x_plane <= p%x_Full ) then + ScBuff_loc = 1.0_ReKi + else if ( p%d_Buff > 0.0_ReKi ) then + ScBuff_loc = max( ( p%x_Buff - x_plane ) / p%d_Buff , 0.0_ReKi ) + else + ScBuff_loc = 0.0_ReKi + end if + end function BufferScale + subroutine Calc_k_WAT() integer(intKi) :: i, iy, iz real(ReKi) :: C, S, dvdr, dvdtheta_r, R, r_tmp @@ -1850,47 +1922,58 @@ subroutine InitStatesWithInputs(numPlanes, numRadii, u, p, xd, m, errStat, errMs character(ErrMsgLen) :: ErrMsg2 real(ReKi) :: correction(3) ! Note, all of these states will have been set to zero in the WD_Init routine - - + ErrStat = ErrID_None ErrMsg = "" - - + correction = 0.0_ReKi do i = 0, 1 xd%x_plane (i) = u%Vx_rel_disk*real(i,ReKi)*real(p%DT_low,ReKi) xd%YawErr_filt (i) = u%YawErr xd%psi_skew_filt = u%psi_skew xd%chi_skew_filt = u%chi_skew - + correction = correction + GetYawCorrection(u%YawErr, u%xhat_disk, xd%x_plane(i), p, errStat2, errMsg2) call SetErrStat(errStat2, errMsg2, errStat, errMsg, RoutineName) if (errStat >= AbortErrLev) then ! TEST: E3 return end if - + !correction = ( p%C_HWkDfl_x + p%C_HWkDfl_xY*u%YawErr )*xd%x_plane(i) + correctionA - + xd%p_plane (:,i) = u%p_hub(:) + xd%x_plane(i)*u%xhat_disk(:) + correction xd%xhat_plane(:,i) = u%xhat_disk(:) xd%V_plane_filt(:,i) = u%V_plane(:,i) xd%Vx_wind_disk_filt(i) = u%Vx_wind_disk xd%TI_amb_filt (i) = u%TI_amb xd%D_rotor_filt (i) = u%D_rotor - - end do - + xd%Vx_rel_disk_filt = u%Vx_rel_disk - + ! Initialze Ct_azavg_filt, Cq_azavg_filt, and Vx_wake; Vr_wake is already initialized to zero, so, we don't need to do that here. xd%Ct_azavg_filt (:) = u%Ct_azavg(:) xd%Cq_azavg_filt (:) = u%Cq_azavg(:) - + call NearWakeCorrection( xd%Ct_azavg_filt, xd%Cq_azavg_filt, xd%Vx_rel_disk_filt, p, m, xd%Vx_wake(:,0), m%Vt_wake, xd%D_rotor_filt(0), errStat, errMsg ) xd%Vx_wake(:,1) = xd%Vx_wake(:,0) - + + + ! Initialize states for cartesian and curled wake formulations + if (p%Mod_Wake == Mod_Wake_Cartesian .or. p%Mod_Wake == Mod_Wake_Curl) then + ! Compute Vx(r) + call NearWakeCorrection( xd%Ct_azavg_filt, xd%Cq_azavg_filt, xd%Vx_rel_disk_filt, p, m, m%Vx_polar(:), m%Vt_wake, xd%D_rotor_filt(0), errStat2, errMsg2 ) + call SetErrStat(ErrStat2, ErrMsg2, errStat, errMsg, RoutineName) + if (errStat >= AbortErrLev) return + call Axisymmetric2CartesianVx(m%Vx_polar, p%r, p%y, p%z, xd%Vx_wake2(:,:,0)) + xd%Vy_wake2(:,:,0) = 0.0_ReKi + xd%Vz_wake2(:,:,0) = 0.0_ReKi + xd%Vx_wake2(:,:,1) = xd%Vx_wake2(:,:,0) + xd%Vy_wake2(:,:,1) = xd%Vy_wake2(:,:,0) + xd%Vz_wake2(:,:,1) = xd%Vz_wake2(:,:,0) + endif + end subroutine InitStatesWithInputs !---------------------------------------------------------------------------------------------------------------------------------- diff --git a/modules/wakedynamics/src/WakeDynamics_Registry.txt b/modules/wakedynamics/src/WakeDynamics_Registry.txt index 692d06973b..424ed0cae1 100644 --- a/modules/wakedynamics/src/WakeDynamics_Registry.txt +++ b/modules/wakedynamics/src/WakeDynamics_Registry.txt @@ -76,6 +76,7 @@ typedef ^ InitInputType WD_InputFileType InputFileData - - - "FAST.Farm inpu typedef ^ InitInputType IntKi TurbNum - 0 - "Turbine ID number (start with 1; end with number of turbines)" - typedef ^ InitInputType CHARACTER(1024) OutFileRoot - - - "The root name derived from the primary FAST.Farm input file" - typedef ^ InitInputType IntKi MaxNumPlanes - - - "Maximum number of wake planes allowed" - +typedef ^ InitInputType ReKi LowResBounds {3}{2} - - "Physical bounds of the low-resolution domain; index 1=XYZ, index 2=lower/upper" m # Define outputs from the initialization routine here: typedef ^ InitOutputType CHARACTER(ChanLen) WriteOutputHdr {:} - - "Names of the output-to-file channels" - @@ -113,6 +114,7 @@ typedef ^ ConstraintStateType ReKi DummyConstrState - - - "R # Define any other states, including integer or logical states here: typedef ^ OtherStateType LOGICAL firstPass - - - "Flag indicating whether or not the states have been initialized with proper inputs" - +typedef ^ OtherStateType LOGICAL MaxPlanesWarned - - - "Flag indicating the MaxNumPlanes-exceeded warning has already been issued for this turbine" - # ..... Misc/Optimization variables................................................................................................. # Define any data that are used only for efficiency purposes (these variables are not associated with time): @@ -186,6 +188,7 @@ typedef ^ ParameterType Logical OutAllPlanes - - - "Output all pl typedef ^ ParameterType CHARACTER(1024) OutFileRoot - - - "The root name derived from the primary FAST.Farm input file" - typedef ^ ParameterType CHARACTER(1024) OutFileVTKDir - - - "The parent directory for all VTK files written by WD" - typedef ^ ParameterType IntKi TurbNum - 0 - "Turbine ID number (start with 1; end with number of turbines)" - +typedef ^ ParameterType ReKi LowResBounds {3}{2} - - "Physical bounds of the low-resolution domain; index 1=XYZ, index 2=lower/upper" m # wake added turbulence (WAT) parameters typedef ^ ParameterType Logical WAT - - - "Switch for turning on and off wake-added turbulence" - typedef ^ ParameterType ReKi WAT_k_Def_k_c - - - "Calibrated parameter for the influence of the maximum wake deficit on wake-added turblence (-) [>=0] or DEFAULT [DEFAULT=0.6]" - diff --git a/modules/wakedynamics/src/WakeDynamics_Types.f90 b/modules/wakedynamics/src/WakeDynamics_Types.f90 index faec69b17a..a1a37b4403 100644 --- a/modules/wakedynamics/src/WakeDynamics_Types.f90 +++ b/modules/wakedynamics/src/WakeDynamics_Types.f90 @@ -92,6 +92,7 @@ MODULE WakeDynamics_Types INTEGER(IntKi) :: TurbNum = 0 !< Turbine ID number (start with 1; end with number of turbines) [-] CHARACTER(1024) :: OutFileRoot !< The root name derived from the primary FAST.Farm input file [-] INTEGER(IntKi) :: MaxNumPlanes = 0_IntKi !< Maximum number of wake planes allowed [-] + REAL(ReKi) , DIMENSION(1:3,1:2) :: LowResBounds = 0.0_ReKi !< Physical bounds of the low-resolution domain; index 1=XYZ, index 2=lower/upper [m] END TYPE WD_InitInputType ! ======================= ! ========= WD_InitOutputType ======= @@ -137,6 +138,7 @@ MODULE WakeDynamics_Types ! ========= WD_OtherStateType ======= TYPE, PUBLIC :: WD_OtherStateType LOGICAL :: firstPass = .false. !< Flag indicating whether or not the states have been initialized with proper inputs [-] + LOGICAL :: MaxPlanesWarned = .false. !< Flag indicating the MaxNumPlanes-exceeded warning has already been issued for this turbine [-] END TYPE WD_OtherStateType ! ======================= ! ========= WD_MiscVarType ======= @@ -208,6 +210,7 @@ MODULE WakeDynamics_Types CHARACTER(1024) :: OutFileRoot !< The root name derived from the primary FAST.Farm input file [-] CHARACTER(1024) :: OutFileVTKDir !< The parent directory for all VTK files written by WD [-] INTEGER(IntKi) :: TurbNum = 0 !< Turbine ID number (start with 1; end with number of turbines) [-] + REAL(ReKi) , DIMENSION(1:3,1:2) :: LowResBounds = 0.0_ReKi !< Physical bounds of the low-resolution domain; index 1=XYZ, index 2=lower/upper [m] LOGICAL :: WAT = .false. !< Switch for turning on and off wake-added turbulence [-] REAL(ReKi) :: WAT_k_Def_k_c = 0.0_ReKi !< Calibrated parameter for the influence of the maximum wake deficit on wake-added turblence (-) [>=0] or DEFAULT [DEFAULT=0.6] [-] REAL(ReKi) :: WAT_k_Def_FMin = 0.0_ReKi !< Calibrated parameter in the eddy viscosity filter function for the WAT maximum wake deficit defining the value in the minimum region [>=0.0 and <=1.0] or DEFAULT [DEFAULT=0.0] [-] @@ -457,6 +460,7 @@ subroutine WD_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode, ErrSta DstInitInputData%TurbNum = SrcInitInputData%TurbNum DstInitInputData%OutFileRoot = SrcInitInputData%OutFileRoot DstInitInputData%MaxNumPlanes = SrcInitInputData%MaxNumPlanes + DstInitInputData%LowResBounds = SrcInitInputData%LowResBounds end subroutine subroutine WD_DestroyInitInput(InitInputData, ErrStat, ErrMsg) @@ -481,6 +485,7 @@ subroutine WD_PackInitInput(RF, Indata) call RegPack(RF, InData%TurbNum) call RegPack(RF, InData%OutFileRoot) call RegPack(RF, InData%MaxNumPlanes) + call RegPack(RF, InData%LowResBounds) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -493,6 +498,7 @@ subroutine WD_UnPackInitInput(RF, OutData) call RegUnpack(RF, OutData%TurbNum); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%OutFileRoot); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%MaxNumPlanes); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%LowResBounds); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine WD_CopyInitOutput(SrcInitOutputData, DstInitOutputData, CtrlCode, ErrStat, ErrMsg) @@ -972,6 +978,7 @@ subroutine WD_CopyOtherState(SrcOtherStateData, DstOtherStateData, CtrlCode, Err ErrStat = ErrID_None ErrMsg = '' DstOtherStateData%firstPass = SrcOtherStateData%firstPass + DstOtherStateData%MaxPlanesWarned = SrcOtherStateData%MaxPlanesWarned end subroutine subroutine WD_DestroyOtherState(OtherStateData, ErrStat, ErrMsg) @@ -989,6 +996,7 @@ subroutine WD_PackOtherState(RF, Indata) character(*), parameter :: RoutineName = 'WD_PackOtherState' if (RF%ErrStat >= AbortErrLev) return call RegPack(RF, InData%firstPass) + call RegPack(RF, InData%MaxPlanesWarned) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -998,6 +1006,7 @@ subroutine WD_UnPackOtherState(RF, OutData) character(*), parameter :: RoutineName = 'WD_UnPackOtherState' if (RF%ErrStat /= ErrID_None) return call RegUnpack(RF, OutData%firstPass); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%MaxPlanesWarned); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine WD_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) @@ -1409,7 +1418,7 @@ subroutine WD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) integer(IntKi), intent(in ) :: CtrlCode integer(IntKi), intent( out) :: ErrStat character(*), intent( out) :: ErrMsg - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: LB(2), UB(2) integer(IntKi) :: ErrStat2 character(*), parameter :: RoutineName = 'WD_CopyParam' ErrStat = ErrID_None @@ -1487,6 +1496,7 @@ subroutine WD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%OutFileRoot = SrcParamData%OutFileRoot DstParamData%OutFileVTKDir = SrcParamData%OutFileVTKDir DstParamData%TurbNum = SrcParamData%TurbNum + DstParamData%LowResBounds = SrcParamData%LowResBounds DstParamData%WAT = SrcParamData%WAT DstParamData%WAT_k_Def_k_c = SrcParamData%WAT_k_Def_k_c DstParamData%WAT_k_Def_FMin = SrcParamData%WAT_k_Def_FMin @@ -1563,6 +1573,7 @@ subroutine WD_PackParam(RF, Indata) call RegPack(RF, InData%OutFileRoot) call RegPack(RF, InData%OutFileVTKDir) call RegPack(RF, InData%TurbNum) + call RegPack(RF, InData%LowResBounds) call RegPack(RF, InData%WAT) call RegPack(RF, InData%WAT_k_Def_k_c) call RegPack(RF, InData%WAT_k_Def_FMin) @@ -1581,7 +1592,7 @@ subroutine WD_UnPackParam(RF, OutData) type(RegFile), intent(inout) :: RF type(WD_ParameterType), intent(inout) :: OutData character(*), parameter :: RoutineName = 'WD_UnPackParam' - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: LB(2), UB(2) integer(IntKi) :: stat logical :: IsAllocAssoc if (RF%ErrStat /= ErrID_None) return @@ -1625,6 +1636,7 @@ subroutine WD_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%OutFileRoot); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%OutFileVTKDir); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%TurbNum); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%LowResBounds); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WAT); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WAT_k_Def_k_c); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WAT_k_Def_FMin); if (RegCheckErr(RF, RoutineName)) return diff --git a/modules/wakedynamics/tests/test_addvelocitycurl.F90 b/modules/wakedynamics/tests/test_addvelocitycurl.F90 new file mode 100644 index 0000000000..954ea5d04d --- /dev/null +++ b/modules/wakedynamics/tests/test_addvelocitycurl.F90 @@ -0,0 +1,42 @@ +module test_addvelocitycurl + +use testdrive, only: new_unittest, unittest_type, error_type, check +use WakeDynamics +use NWTC_Library + +implicit none +private +public :: test_addvelocitycurl_suite + +contains + +!> Collect all exported unit tests +subroutine test_addvelocitycurl_suite(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + testsuite = [ & + new_unittest("test_add_velocity_curl", test_add_velocity_curl) & + ] +end subroutine + +!> Checks the curled-wake velocity-curl calculation against known reference values. +subroutine test_add_velocity_curl(error) + type(error_type), allocatable, intent(out) :: error + + real(ReKi) :: Vy_curl(2, 2) = 0.0_ReKi + real(ReKi) :: Vz_curl(2, 2) = 0.0_ReKi + real(ReKi) :: y(2) = (/0., 2./) + real(ReKi) :: z(2) = (/-1., 1./) + real(ReKi) :: Gamma0 + + call AddVelocityCurl(Vx=10., yaw_angle=0.1, nVortex=100, R=63., psi_skew=0.2, & + y=y, z=z, Ct_avg=0.7, sigma_d=0.2, Vy_curl=Vy_curl, Vz_curl=Vz_curl, Gamma0=Gamma0) + + call check(error, abs(Vy_curl(1, 1) + 0.217109) < 1e-4, "Vy_curl(1,1) does not match reference value") + if (allocated(error)) return + + call check(error, abs(Vz_curl(2, 2) + 4.459746e-2) < 1e-4, "Vz_curl(2,2) does not match reference value") + if (allocated(error)) return + +end subroutine + +end module diff --git a/modules/wakedynamics/tests/test_axisymmetric2cartesian.F90 b/modules/wakedynamics/tests/test_axisymmetric2cartesian.F90 new file mode 100644 index 0000000000..a4bbc509ed --- /dev/null +++ b/modules/wakedynamics/tests/test_axisymmetric2cartesian.F90 @@ -0,0 +1,58 @@ +module test_axisymmetric2cartesian + +use testdrive, only: new_unittest, unittest_type, error_type, check +use WakeDynamics +use NWTC_Library + +implicit none +private +public :: test_axisymmetric2cartesian_suite + +contains + +!> Collect all exported unit tests +subroutine test_axisymmetric2cartesian_suite(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + testsuite = [ & + new_unittest("test_axisymmetric2cartesian_vel", test_axisymmetric2cartesian_vel) & + ] +end subroutine + +!> Checks that converting an axisymmetric radial velocity field to Cartesian coordinates +!! recovers the known radial velocity magnitude and axial velocity at each grid point. +subroutine test_axisymmetric2cartesian_vel(error) + type(error_type), allocatable, intent(out) :: error + + real(ReKi) :: r(4) = (/0., 1., 2., 3./) + real(ReKi) :: y(4) = (/0., 1., 1.5, 2./) + real(ReKi) :: z(5) = (/0., 0.5, 1., 1.5, 2./) + real(ReKi) :: Vr_axi(4) + real(ReKi) :: Vx_axi(4) + real(ReKi) :: Vx(4, 5) = 0.0_ReKi + real(ReKi) :: Vy(4, 5) = 0.0_ReKi + real(ReKi) :: Vz(4, 5) = 0.0_ReKi + integer :: i, j + real(ReKi) :: Vr, r_tmp + character(100) :: label + + Vr_axi = 4._ReKi*r + Vx_axi = 3._ReKi*r + call Axisymmetric2CartesianVel(Vx_axi, Vr_axi, r, y, z, Vx, Vy, Vz) + + do i = 1, size(y) + do j = 1, size(z) + r_tmp = sqrt(y(i)**2 + z(j)**2) + Vr = sqrt(Vy(i, j)**2 + Vz(i, j)**2) + + write (label, '(A,I0,A,I0,A)') "Vr mismatch at (", i, ",", j, ")" + call check(error, abs(Vr - 4*r_tmp) < 1e-3, trim(label)) + if (allocated(error)) return + + write (label, '(A,I0,A,I0,A)') "Vx mismatch at (", i, ",", j, ")" + call check(error, abs(Vx(i, j) - 3*r_tmp) < 1e-3, trim(label)) + if (allocated(error)) return + end do + end do +end subroutine + +end module diff --git a/modules/wakedynamics/tests/test_wd_maxplanes_warning.F90 b/modules/wakedynamics/tests/test_wd_maxplanes_warning.F90 new file mode 100644 index 0000000000..3274580582 --- /dev/null +++ b/modules/wakedynamics/tests/test_wd_maxplanes_warning.F90 @@ -0,0 +1,117 @@ +module test_wd_maxplanes_warning + +use testdrive, only: new_unittest, unittest_type, error_type, check +use WakeDynamics +use WakeDynamics_Types +use NWTC_Library + +implicit none +private +public :: test_wd_maxplanes_warning_suite + +contains + +!> Collect all exported unit tests +subroutine test_wd_maxplanes_warning_suite(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + testsuite = [ & + new_unittest("test_warning_issued_only_once", test_warning_issued_only_once) & + ] +end subroutine + +!> Regression test: once the number of tracked wake planes exceeds p%MaxNumPlanes, +!! WD_UpdateStates clamps the count and issues an ErrID_Warn. With a small +!! MaxNumPlanes and no removal (merging/buffer-exit) happening, this condition +!! recurs on every subsequent call. OtherState%MaxPlanesWarned should ensure the +!! warning is only issued on the first occurrence, not on every following call. +subroutine test_warning_issued_only_once(error) + type(error_type), allocatable, intent(out) :: error + + type(WD_InitInputType) :: InitInp + type(WD_InputType) :: u + type(WD_ParameterType) :: p + type(WD_ContinuousStateType) :: x + type(WD_DiscreteStateType) :: xd + type(WD_ConstraintStateType) :: z + type(WD_OtherStateType) :: OtherState + type(WD_OutputType) :: y + type(WD_MiscVarType) :: m + type(WD_InitOutputType) :: InitOut + integer(IntKi) :: errStat + character(ErrMsgLen) :: errMsg + real(DbKi), parameter :: DT_low = 1.0_DbKi + integer(IntKi), parameter :: MaxNumPlanes = 3 + integer(IntKi) :: n + integer(IntKi) :: warnCount + + errStat = ErrID_None + errMsg = '' + + ! --- Minimal, valid WakeDynamics (Polar model) initialization input + InitInp%TurbNum = 1 + InitInp%OutFileRoot = 'test_wd_maxplanes_warning' + InitInp%MaxNumPlanes = MaxNumPlanes + ! Generous low-res domain bounds: planes are never dropped via OOB merging, + ! so the only thing keeping NumPlanes at the cap is the clamp itself. + InitInp%LowResBounds(:, 1) = -1.0e6_ReKi + InitInp%LowResBounds(:, 2) = 1.0e6_ReKi + + InitInp%InputFileData%Mod_Wake = Mod_Wake_Polar + InitInp%InputFileData%RotorDiamRef = 20.0_ReKi + InitInp%InputFileData%dr = 5.0_ReKi + InitInp%InputFileData%NumRadii = 3 + InitInp%InputFileData%NumDFull = 15.0_ReKi + InitInp%InputFileData%NumDBuff = 5.0_ReKi + InitInp%InputFileData%f_c = 0.17_ReKi + InitInp%InputFileData%C_NearWake = 1.8_ReKi + InitInp%InputFileData%k_vAmb = 0.05_ReKi + InitInp%InputFileData%C_vAmb_FMin = 1.0_ReKi + InitInp%InputFileData%C_vAmb_DMin = 0.0_ReKi + InitInp%InputFileData%C_vAmb_DMax = 1.0_ReKi + InitInp%InputFileData%C_vAmb_Exp = 0.01_ReKi + InitInp%InputFileData%k_vShr = 0.016_ReKi + InitInp%InputFileData%C_vShr_FMin = 0.2_ReKi + InitInp%InputFileData%C_vShr_DMin = 3.0_ReKi + InitInp%InputFileData%C_vShr_DMax = 25.0_ReKi + InitInp%InputFileData%C_vShr_Exp = 0.1_ReKi + InitInp%InputFileData%Mod_WakeDiam = 1 + + call WD_Init(InitInp, u, p, x, xd, z, OtherState, y, m, DT_low, InitOut, errStat, errMsg) + call check(error, errStat < AbortErrLev, "WD_Init failed: "//trim(errMsg)) + if (allocated(error)) return + + ! --- Steady, uniform axial inflow so plane advection is simple and predictable + u%xhat_disk = (/1.0_ReKi, 0.0_ReKi, 0.0_ReKi/) + u%YawErr = 0.0_ReKi + u%psi_skew = 0.0_ReKi + u%chi_skew = 0.0_ReKi + u%p_hub = (/0.0_ReKi, 0.0_ReKi, 90.0_ReKi/) + u%Vx_wind_disk = 8.0_ReKi + u%TI_amb = 0.1_ReKi + u%D_rotor = 20.0_ReKi + u%Vx_rel_disk = 8.0_ReKi + u%Ct_azavg = 0.0_ReKi + u%Cq_azavg = 0.0_ReKi + u%V_plane(1, :) = 8.0_ReKi + u%V_plane(2, :) = 0.0_ReKi + u%V_plane(3, :) = 0.0_ReKi + + ! --- Advance well past MaxNumPlanes so the exceeded-cap clamp is hit repeatedly. + warnCount = 0 + do n = 0, 4*MaxNumPlanes + errStat = ErrID_None + errMsg = '' + call WD_UpdateStates(real(n, DbKi)*DT_low, n, u, p, x, xd, z, OtherState, m, errStat, errMsg) + call check(error, errStat < AbortErrLev, "WD_UpdateStates failed: "//trim(errMsg)) + if (allocated(error)) return + if (index(errMsg, 'exceeded the allowed number') > 0) warnCount = warnCount + 1 + call check(error, NINT(xd%NumPlanes) <= p%MaxNumPlanes, "NumPlanes exceeded MaxNumPlanes despite clamp") + if (allocated(error)) return + end do + + call check(error, warnCount == 1, "Expected the MaxNumPlanes-exceeded warning exactly once") + if (allocated(error)) return + +end subroutine + +end module diff --git a/modules/wakedynamics/tests/test_wd_oobidx.F90 b/modules/wakedynamics/tests/test_wd_oobidx.F90 new file mode 100644 index 0000000000..2d59a3192d --- /dev/null +++ b/modules/wakedynamics/tests/test_wd_oobidx.F90 @@ -0,0 +1,124 @@ +module test_wd_oobidx + +use testdrive, only: new_unittest, unittest_type, error_type, check +use WakeDynamics +use WakeDynamics_Types +use NWTC_Library + +implicit none +private +public :: test_wd_oobidx_suite + +contains + +!> Collect all exported unit tests +subroutine test_wd_oobidx_suite(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + testsuite = [ & + new_unittest("test_all_planes_simultaneously_out_of_bounds", test_all_planes_simultaneously_out_of_bounds) & + ] +end subroutine + +!> Regression test for a stack array-bounds overflow in WD_UpdateStates: when every +!! existing wake plane (up to p%MaxNumPlanes of them) is simultaneously found to be +!! outside the low-resolution domain bounds, the local automatic array `oobIdx` must +!! have room to record all of them. Before the fix, `oobIdx` was declared +!! `0:MaxNumPlanes-1` (size MaxNumPlanes) -- one element too small to hold +!! `MaxNumPlanes` out-of-bounds indices -- causing an out-of-bounds stack write. +!! (With gfortran's `-fcheck=all` Debug flag, this manifests as a hard runtime +!! crash rather than a graceful test failure.) +subroutine test_all_planes_simultaneously_out_of_bounds(error) + type(error_type), allocatable, intent(out) :: error + + type(WD_InitInputType) :: InitInp + type(WD_InputType) :: u + type(WD_ParameterType) :: p + type(WD_ContinuousStateType) :: x + type(WD_DiscreteStateType) :: xd + type(WD_ConstraintStateType) :: z + type(WD_OtherStateType) :: OtherState + type(WD_OutputType) :: y + type(WD_MiscVarType) :: m + type(WD_InitOutputType) :: InitOut + integer(IntKi) :: errStat + character(ErrMsgLen) :: errMsg + real(DbKi), parameter :: DT_low = 1.0_DbKi + integer(IntKi), parameter :: MaxNumPlanes = 4 + integer(IntKi) :: n + + errStat = ErrID_None + errMsg = '' + + ! --- Minimal, valid WakeDynamics (Polar model) initialization input + InitInp%TurbNum = 1 + InitInp%OutFileRoot = 'test_wd_oobidx' + InitInp%MaxNumPlanes = MaxNumPlanes + ! Generous low-res domain bounds so no plane is dropped during natural growth + InitInp%LowResBounds(:, 1) = -1.0e6_ReKi + InitInp%LowResBounds(:, 2) = 1.0e6_ReKi + + InitInp%InputFileData%Mod_Wake = Mod_Wake_Polar + InitInp%InputFileData%RotorDiamRef = 20.0_ReKi + InitInp%InputFileData%dr = 5.0_ReKi + InitInp%InputFileData%NumRadii = 3 + InitInp%InputFileData%NumDFull = 15.0_ReKi + InitInp%InputFileData%NumDBuff = 5.0_ReKi + InitInp%InputFileData%f_c = 0.17_ReKi + InitInp%InputFileData%C_NearWake = 1.8_ReKi + InitInp%InputFileData%k_vAmb = 0.05_ReKi + InitInp%InputFileData%C_vAmb_FMin = 1.0_ReKi + InitInp%InputFileData%C_vAmb_DMin = 0.0_ReKi + InitInp%InputFileData%C_vAmb_DMax = 1.0_ReKi + InitInp%InputFileData%C_vAmb_Exp = 0.01_ReKi + InitInp%InputFileData%k_vShr = 0.016_ReKi + InitInp%InputFileData%C_vShr_FMin = 0.2_ReKi + InitInp%InputFileData%C_vShr_DMin = 3.0_ReKi + InitInp%InputFileData%C_vShr_DMax = 25.0_ReKi + InitInp%InputFileData%C_vShr_Exp = 0.1_ReKi + InitInp%InputFileData%Mod_WakeDiam = 1 + + call WD_Init(InitInp, u, p, x, xd, z, OtherState, y, m, DT_low, InitOut, errStat, errMsg) + call check(error, errStat < AbortErrLev, "WD_Init failed: "//trim(errMsg)) + if (allocated(error)) return + + ! --- Steady, uniform axial inflow so plane advection is simple and predictable + u%xhat_disk = (/1.0_ReKi, 0.0_ReKi, 0.0_ReKi/) + u%YawErr = 0.0_ReKi + u%psi_skew = 0.0_ReKi + u%chi_skew = 0.0_ReKi + u%p_hub = (/0.0_ReKi, 0.0_ReKi, 90.0_ReKi/) + u%Vx_wind_disk = 8.0_ReKi + u%TI_amb = 0.1_ReKi + u%D_rotor = 20.0_ReKi + u%Vx_rel_disk = 8.0_ReKi + u%Ct_azavg = 0.0_ReKi + u%Cq_azavg = 0.0_ReKi + u%V_plane(1, :) = 8.0_ReKi + u%V_plane(2, :) = 0.0_ReKi + u%V_plane(3, :) = 0.0_ReKi + + ! --- Grow the number of tracked wake planes up to p%MaxNumPlanes + do n = 0, MaxNumPlanes - 3 + call WD_UpdateStates(real(n, DbKi)*DT_low, n, u, p, x, xd, z, OtherState, m, errStat, errMsg) + call check(error, errStat < AbortErrLev, "WD_UpdateStates failed during plane growth: "//trim(errMsg)) + if (allocated(error)) return + end do + call check(error, NINT(xd%NumPlanes) == p%MaxNumPlanes, "Expected wake planes to grow to MaxNumPlanes") + if (allocated(error)) return + + ! --- Tighten the low-resolution domain bounds so every existing plane + ! (including the disk plane, which stays near Z=90) is now out of bounds. + p%LowResBounds(:, 1) = 500.0_ReKi + p%LowResBounds(:, 2) = 600.0_ReKi + + ! With the fix, this call completes without an out-of-bounds write to the + ! oobIdx local array (previously sized one element too small to hold + ! MaxNumPlanes out-of-bounds indices). + call WD_UpdateStates(real(MaxNumPlanes - 2, DbKi)*DT_low, MaxNumPlanes - 2, u, p, x, xd, z, OtherState, m, errStat, errMsg) + call check(error, errStat < AbortErrLev, & + "WD_UpdateStates failed when all planes were simultaneously out of bounds: "//trim(errMsg)) + if (allocated(error)) return + +end subroutine + +end module diff --git a/modules/wakedynamics/tests/wakedynamics_utest.F90 b/modules/wakedynamics/tests/wakedynamics_utest.F90 new file mode 100644 index 0000000000..d131437d9c --- /dev/null +++ b/modules/wakedynamics/tests/wakedynamics_utest.F90 @@ -0,0 +1,39 @@ +program wakedynamics_utest +use, intrinsic :: iso_fortran_env, only: error_unit +use testdrive, only: run_testsuite, new_testsuite, testsuite_type + +use test_addvelocitycurl, only: test_addvelocitycurl_suite +use test_axisymmetric2cartesian, only: test_axisymmetric2cartesian_suite +use test_wd_oobidx, only: test_wd_oobidx_suite +use test_wd_maxplanes_warning, only: test_wd_maxplanes_warning_suite +use NWTC_Num + +implicit none +integer :: stat, is +type(testsuite_type), allocatable :: testsuites(:) +character(len=*), parameter :: fmt = '("#", *(1x, a))' + +stat = 0 + +call SetConstants() + +testsuites = [ & + new_testsuite("AddVelocityCurl", test_addvelocitycurl_suite), & + new_testsuite("Axisymmetric2Cartesian", test_axisymmetric2cartesian_suite), & + new_testsuite("WD_UpdateStates_OOBIdx", test_wd_oobidx_suite), & + new_testsuite("WD_UpdateStates_MaxPlanesWarning", test_wd_maxplanes_warning_suite) & + ] + +do is = 1, size(testsuites) + write (error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat, parallel=.false.) +end do + +if (stat > 0) then + write (error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop +end if + +write (error_unit, fmt) "All tests PASSED" + +end program diff --git a/reg_tests/r-test b/reg_tests/r-test index 2aa1664832..861afeb59b 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 2aa1664832f47d6c1286583b33c36bf4dfe1b80b +Subproject commit 861afeb59b3b7b9c59b2f205ef2faf65bcb25157 diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt index c8847a0957..52f4e1c7f1 100644 --- a/unit_tests/CMakeLists.txt +++ b/unit_tests/CMakeLists.txt @@ -62,6 +62,17 @@ add_executable(inflowwind_utest target_link_libraries(inflowwind_utest ifwlib versioninfolib testdrivelib) add_test(NAME inflowwind_utest COMMAND inflowwind_utest) +# WakeDynamics Unit Tests +add_executable(wakedynamics_utest + ${PROJECT_SOURCE_DIR}/modules/wakedynamics/tests/wakedynamics_utest.F90 + ${PROJECT_SOURCE_DIR}/modules/wakedynamics/tests/test_addvelocitycurl.F90 + ${PROJECT_SOURCE_DIR}/modules/wakedynamics/tests/test_axisymmetric2cartesian.F90 + ${PROJECT_SOURCE_DIR}/modules/wakedynamics/tests/test_wd_oobidx.F90 + ${PROJECT_SOURCE_DIR}/modules/wakedynamics/tests/test_wd_maxplanes_warning.F90 +) +target_link_libraries(wakedynamics_utest wdlib versioninfolib testdrivelib) +add_test(NAME wakedynamics_utest COMMAND wakedynamics_utest) + # NWTC Library Unit Tests add_executable(nwtc_library_utest ${PROJECT_SOURCE_DIR}/modules/nwtc-library/tests/nwtc_library_utest.F90 @@ -92,6 +103,7 @@ add_custom_target( beamdyn_utest aerodyn_utest inflowwind_utest + wakedynamics_utest ) if (AMREX_READER AND BUILD_FASTFARM) @@ -105,6 +117,10 @@ if (AMREX_READER AND BUILD_FASTFARM) ${CMAKE_SOURCE_DIR}/modules/awae/tests/test_AMReX_reader.f90 ) target_link_libraries(amrex_reader_utest awaelib nwtclibs testdrivelib) + # Both sources are Fortran and main is a Fortran program, but awaelib pulls in the C++ + # awaelib_c, so CMake picks the C++ driver to link with and its crt1.o then looks for a + # C main. Pin the linker to Fortran, as glue-codes/fast-farm does for the same reason. + set_target_properties(amrex_reader_utest PROPERTIES LINKER_LANGUAGE Fortran) add_test(NAME amrex_reader_utest COMMAND amrex_reader_utest) add_dependencies(amrex_reader_utest copy_amrex_reader_utest_files) # copy input files needed for this test add_dependencies(unit_tests amrex_reader_utest) # add to the unit_tests target diff --git a/vs-build/modules/AWAE.vfproj b/vs-build/modules/AWAE.vfproj index c162009a2e..35dcc383d3 100644 --- a/vs-build/modules/AWAE.vfproj +++ b/vs-build/modules/AWAE.vfproj @@ -95,6 +95,7 @@ +