diff --git a/.gitignore b/.gitignore index 09166cd156..74be88339b 100644 --- a/.gitignore +++ b/.gitignore @@ -64,3 +64,8 @@ varcache # Python cache files openfast_io/dist/ openfast_io/openfast_io/_version.py + +# docker with VScode config +*.code-workspace +docker-compose.yml +.devcontainer/ 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/api_change.rst b/docs/source/user/api_change.rst index 9a2087ff80..8db1e466f5 100644 --- a/docs/source/user/api_change.rst +++ b/docs/source/user/api_change.rst @@ -9,6 +9,31 @@ The changes are tabulated according to the module input file, line number, and f The line number corresponds to the resulting line number after all changes are implemented. Thus, be sure to implement each in order so that subsequent line numbers are correct. +OpenFAST v5.0.x to next release +------------------------------- + +Two new optional VISUALIZATION sub-sections were added to the FAST.Farm +primary input file to give users finer control over what gets written +as VTK output. Legacy decks that do not include either block continue +to work unchanged. + +The first block, ``--- AXIS-ALIGNED PLANE SLICES (extent-controlled) ---``, +lets the user specify axis-aligned uniform-grid slices with explicit +2-D extents anywhere in the low-resolution domain (Feature 2, plan +§3). Output is XML VTK ``.vts`` + ``.vts.series``. The block header +is optional and, when present, sits between the existing VISUALIZATION +block and the OUTPUT block. + +The second block, ``--- TERRAIN-FOLLOWING SAMPLING ---``, lets the user +emit point clouds sampled from STL surfaces or plain-text/CSV files, +lifted to one or more offset sheets along either an explicit normal +or the STL per-facet normal (Feature 3, plan §2). Output is XML VTK +``.vtp`` + ``.vtp.series``. + +See :numref:`FF:Input:PlaneSlices` and :numref:`FF:Input:TerrainSlices` +for the full grammar, and :numref:`FF:sec:SliceOutputs` for guidance +on which slice path to pick. + OpenFAST v4.2.x to OpenFAST v5.0.0 ----------------------------------- diff --git a/docs/source/user/fast.farm/FutureWork.rst b/docs/source/user/fast.farm/FutureWork.rst index 34ecdc3179..09508e5e74 100644 --- a/docs/source/user/fast.farm/FutureWork.rst +++ b/docs/source/user/fast.farm/FutureWork.rst @@ -116,3 +116,12 @@ releases: - Adopt the capability to support undersea marine turbine arrays (which may require supporting direct feedthrough of input to output to handle the added-mass effects). + +- Extend the terrain-following sampling + (:numref:`FF:Input:TerrainSlices`) to also accept arbitrary + non-axis-aligned planar samples. Today the extent-controlled plane + slice feature (:numref:`FF:Input:PlaneSlices`) only accepts the + keywords ``XY``, ``YZ``, and ``XZ`` with a + fatal init error for any other value. See GitHub + issue `#2383 `__ + for the motivating discussion and one proposed grammar. diff --git a/docs/source/user/fast.farm/InputFiles.rst b/docs/source/user/fast.farm/InputFiles.rst index bcc524f4e7..fbff5e8c49 100644 --- a/docs/source/user/fast.farm/InputFiles.rst +++ b/docs/source/user/fast.farm/InputFiles.rst @@ -846,6 +846,182 @@ Therefore, disabling visualization is recommended when running many FAST.Farm simulations. See :numref:`FF:Output:Vis` for visualization output file details. +.. _FF:Input:PlaneSlices: + +Axis-Aligned Plane Slices (extent-controlled) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This optional block lets the user emit an arbitrary number of +axis-aligned uniform-grid VTK slices, each with a user-controlled 2D +extent placed anywhere in the low-resolution domain. Unlike the classic +**NOutDisWindXY** / **YZ** / **XZ** slices above (which always span the +entire low-res domain), each slice here is bounded by an explicit +2-D extent, so a single hub-height cut through a large farm can be +kept to just the region of interest. + +The block is optional. Legacy decks omit it entirely; older files still +parse without modification. + +The block begins with the section header +``--- AXIS-ALIGNED PLANE SLICES (extent-controlled) ---`` and is +followed by: + +- **NumPlaneSlices** [integer] specifies the number of + extent-controlled axis-aligned planar slices to output + (:math:`0` to :math:`99`). Set to :math:`0` to disable the feature. + +- **WrPlaneDT** [sec] specifies the sampling period for this + feature. **DEFAULT** falls back to **DT_Low**. Internally rounded + to the nearest multiple of **DT_Low**. Independent of the classic + **WrDisDT** so that (for example) a small hub-height sheet can be + sampled frequently while heavier full-domain outputs sample rarely. + +- A two-line column header followed by **NumPlaneSlices** slice rows + with the following columns: + + **SliceName** — quoted string. Free-form slice name; appears in the + output file names. Example: ``"T1_0D"``. + + **origin(m)** — 3-vector, in metres. Plane corner (not centre) + in the farm-global frame. Example: ``(0 -300 0)``. + + **plane** — keyword. One of ``XY``, ``YZ``, or ``XZ`` + (case-insensitive). Specifies the orientation of the sampling + plane. Any other value is a fatal init error. + + **extent1(m)** — positive scalar, in metres. In-plane extent + along the first axis of the named plane. Example: ``600``. + + **extent2(m)** — positive scalar, in metres. In-plane extent + along the second axis of the named plane. Example: ``400``. + +The (extent1, extent2) axes are ordered by the named plane: +``XY`` :math:`\Rightarrow` (X, Y), +``YZ`` :math:`\Rightarrow` (Y, Z), +``XZ`` :math:`\Rightarrow` (X, Z). + +Grid resolution is fixed to the low-resolution spacing +(``dX_Low``, ``dY_Low``, ``dZ_Low``); no ``npoints`` column is +required. Both parenthesised and bare numeric vectors parse — the +parser strips ``(``, ``)``, and ``,`` before tokenising. + +Output is written as VTK XML ``StructuredGrid`` (``.vts``) plus a +matching ParaView ``.vts.series`` sidecar so time-series playback +works out of the box. Any node whose global-frame position falls +outside the low-resolution domain is written as IEEE quiet NaN; +ParaView automatically masks these. + +Example row (three slices — one XY at hub height, one YZ through a +turbine, one XZ through the farm centre): + +.. code-block:: none + + --- AXIS-ALIGNED PLANE SLICES (extent-controlled) --- + 3 NumPlaneSlices - Number of axis-aligned planar slices (-) + DEFAULT WrPlaneDT - Sampling period (s) or DEFAULT (=DT_Low) + SliceName origin(m) plane extent1(m) extent2(m) + (-) (m,m,m) (-) (m) (m) + "hubXY" (700 700 95) XY 1200 600 + "crossYZ" (1000 700 5) YZ 600 300 + "crossXZ" (700 1000 5) XZ 1200 300 + +Output naming follows the pattern +``.Plane...vts`` with a companion +``.Plane..vts.series``, all under ``vtk_ff/``. + +.. _FF:Input:TerrainSlices: + +Terrain-Following Sampling +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This optional block lets the user emit an arbitrary number of point +clouds sampled from an STL surface or a plain-text/CSV point list. +Each source can be lifted to multiple offset "sheets" along a +user-provided normal, or along the per-facet normal read from an STL +file. Typical use cases: a mountain surface sampled at three +rotor-tip heights, or a hand-authored hub-height sample grid over an +irregular farm layout. + +The block is optional. Legacy decks omit it entirely. + +The block begins with the section header +``--- TERRAIN-FOLLOWING SAMPLING ---`` and is followed by: + +- **NumTerrainSlices** [integer] specifies the number of + terrain-following slices (:math:`0` to :math:`99`). + +- **WrTerrainDT** [sec] specifies the sampling period for this + feature. **DEFAULT** falls back to **DT_Low**. Independent of + both **WrDisDT** and **WrPlaneDT** so that heavy terrain samples + can be emitted less often than lighter axis-aligned slices. + +- A two-line column header followed by **NumTerrainSlices** slice + rows with the following columns: + + **SliceName** — quoted string. Free-form slice name; appears in + output file names. + + **Offsets(m)** — comma- or whitespace-separated list of scalar + displacements (metres) along the **OffsetNormal**. Each value + spawns a separate sheet sharing the source geometry. An empty + list is equivalent to a single ``0``. + + **OffsetNormal** — 3-vector, or the keyword ``default``. Unit + vector along which the offsets are applied. When ``default`` and + **SourceType** is ``STL``, the per-facet normal from the STL is + used instead. An explicit vector is required when **SourceType** + is ``Point``. + + **SourceType** — keyword. ``STL`` reads an ASCII or binary STL + file and uses each facet's three vertices as sample points (no + de-duplication). ``Point`` reads a plain-text or CSV point cloud + file (see the point-cloud file format description below). + + **FileName** — quoted string. Path to the STL or point-cloud + file. Relative paths are resolved against the FAST.Farm primary + input file directory. + +Example row (two slices — a mountain lifted to three offset sheets +and a nine-point hub-height cloud): + +.. code-block:: none + + --- TERRAIN-FOLLOWING SAMPLING --- + 2 NumTerrainSlices - Number of terrain-following slices (-) + DEFAULT WrTerrainDT - Sampling period (s) or DEFAULT (=DT_Low) + SliceName Offsets(m) OffsetNormal SourceType FileName + (-) (m,list) (-,-,-|dflt) (STL|Point) (quoted) + "terr" 50 100 150 default STL "mountain.stl" + "hub" 0.0 (0 0 1) Point "hub_points.txt" + +**Point-cloud file format.** Plain ASCII, one point per line, three +floating-point columns ``x y z`` in metres in the FAST.Farm global +coordinate frame. Delimiter is whitespace or comma (both ``.txt`` and +``.csv`` are handled by the same parser). Comments are lines whose +first non-blank character is ``#`` or ``!``; trailing comments are +also allowed. Blank lines are ignored. An optional header row is +tolerated as long as it fails numeric parsing (this lets ``.csv`` +exports with an ``x,y,z`` header row load unchanged). Column order is +fixed to ``x y z``. Point count is derived from the file — no +``NumPoints`` header is required. Duplicate points are not +de-duplicated (unlike STL vertices). + +Output is written as VTK XML ``PolyData`` (``.vtp``) plus a matching +``.vtp.series`` sidecar. Any sample point falling outside the +low-resolution domain is written as IEEE quiet NaN, and an +``ErrID_Warn`` line is emitted at init reporting the count and +fraction. When a slice's total sample-location count +(``NPts × N_offsets``) is unusually large (currently +:math:`\ge 10\,000\,000`), an ``ErrID_Info`` size advisory is emitted +reporting the estimated RAM cost. No hard cap is imposed — the +NWTC-library allocator already fails cleanly with a byte-count message +if the OS refuses. + +Output naming follows the pattern +``.TerrSlice...vtp`` with a companion +``.TerrSlice..vtp.series``, all under +``vtk_ff/``. + Output ~~~~~~ diff --git a/docs/source/user/fast.farm/ModelGuidance.rst b/docs/source/user/fast.farm/ModelGuidance.rst index 6ac30dd503..f71757c232 100644 --- a/docs/source/user/fast.farm/ModelGuidance.rst +++ b/docs/source/user/fast.farm/ModelGuidance.rst @@ -463,6 +463,53 @@ error universally across all aeroelastic software that use Mann boxes is needed `future work `__. +.. _FF:sec:SliceOutputs: + +Choosing a Visualization Slice Output +------------------------------------- + +FAST.Farm ships three different slice-output paths, each with its own +place. Choose based on how much of the domain you need and how you plan +to load the output in ParaView. + +**Classic slices** (**NOutDisWindXY** / **NOutDisWindYZ** / +**NOutDisWindXZ**, see :numref:`FF:Input:VTK`). These always span the +entire low-res domain and write legacy ASCII ``.vtk`` files with +``STRUCTURED_POINTS`` topology. Cheap to hand-author, but on a large +farm each XY hub-height slice is potentially the full farm cross +section — hundreds of megabytes per step, mostly of empty space. +Recommended: small farms; quick debugging; existing decks that already +use them. + +**Extent-controlled axis-aligned slices** (**NumPlaneSlices**, see +:numref:`FF:Input:PlaneSlices`). Same axis-aligned uniform grid as +the classic slices, but each slice carries its own explicit 2-D extent +so a hub-height sheet can be kept to just the region of interest. +Output is XML ``.vts`` with a ParaView ``.vts.series`` sidecar for +straightforward time-series playback. Recommended for large-farm +runs where you know the region of interest, and for any new deck +that needs axis-aligned slices; use it in place of the classic +**NOutDisWindXY** / **YZ** / **XZ** blocks when disk-footprint or +per-slice extent control matters. + +**Terrain-following point clouds** (**NumTerrainSlices**, see +:numref:`FF:Input:TerrainSlices`). Reads a surface (STL) or a bare +point list from disk and samples along user-provided offset sheets. +This is the only path that supports non-axis-aligned sampling +today — its typical uses are (i) sampling above a real mountain +surface at rotor-tip heights, (ii) sampling at hub height across an +irregular farm footprint, and (iii) validating that arbitrary +points in the domain carry the expected disturbed velocity. Output +is XML ``.vtp`` with a ``.vtp.series`` sidecar. + +All three paths write into the ``vtk_ff/`` subdirectory of the deck +and are throttled independently: **WrDisDT** for the classic path, +**WrPlaneDT** for **NumPlaneSlices**, **WrTerrainDT** for +**NumTerrainSlices**. Each defaults to **DT_Low** if omitted, but +having independent knobs lets you (say) emit a heavy full-domain +slice once per minute while sampling a small hub-height sheet every +second. + .. _FF:sec:DiscRecs: Low- and High-Resolution Domain Discretization diff --git a/docs/source/user/fast.farm/OutputFiles.rst b/docs/source/user/fast.farm/OutputFiles.rst index 4e65e3accf..f9d2bc94be 100644 --- a/docs/source/user/fast.farm/OutputFiles.rst +++ b/docs/source/user/fast.farm/OutputFiles.rst @@ -81,6 +81,55 @@ Each visualization output file follows the same VTK format used for the ambient wind data files for the high-fidelity precursor simulations. See :numref:`FF:AmbWindIfW` for details on the file format. +.. _FF:Output:PlaneSlices: + +Extent-Controlled Axis-Aligned Plane Slices (Feature 2) +------------------------------------------------------- + +If **NumPlaneSlices** is greater than zero in the FAST.Farm primary +input file (see :numref:`FF:Input:PlaneSlices`), one file per slice per +sampled step is written to *vtk_ff/* with naming +*.Plane...vts* where +** is the user-provided free-form name from the deck and +** is a zero-padded step counter. + +Each slice also has a companion ParaView time-series sidecar named +*.Plane..vts.series* — a small JSON file listing +every ``.vts`` step and its simulation time. Loading the ``.series`` +file in ParaView plays the whole slice as an animation. + +The ``.vts`` format is VTK XML *StructuredGrid* with explicit point +coordinates and a 3-component ``Velocity`` field. Nodes whose +global-frame position falls outside the low-resolution domain carry +IEEE quiet NaN in the ``Velocity`` field; ParaView masks these +automatically without any extra filter. + +The sampling period is set by **WrPlaneDT** in the deck (defaults to +**DT_Low**). It is independent of the classic **WrDisDT** so that +lightweight hub-height sheets can be emitted at a much higher rate +than heavier full-domain outputs. + +.. _FF:Output:TerrainSlices: + +Terrain-Following Point Cloud Slices (Feature 3) +------------------------------------------------ + +If **NumTerrainSlices** is greater than zero in the FAST.Farm primary +input file (see :numref:`FF:Input:TerrainSlices`), one file per slice +per sampled step is written to *vtk_ff/* with naming +*.TerrSlice...vtp* and a +companion *.TerrSlice..vtp.series* sidecar. + +The ``.vtp`` format is VTK XML *PolyData* with per-point coordinates +and a 3-component ``Velocity`` field. Each source vertex is +replicated once per offset sheet, so a slice with :math:`N` source +points and :math:`M` offsets writes :math:`N \cdot M` points per +step. Nodes whose sample location falls outside the low-resolution +domain are set to IEEE quiet NaN and masked by ParaView. + +The sampling period is set by **WrTerrainDT** in the deck (defaults to +**DT_Low**). + Visualizing the ambient wind and wake interactions can be useful for interpreting results and debugging problems. However, FAST.Farm will generate many files per output option when **WrDisWind** = TRUE and/or diff --git a/glue-codes/fast-farm/src/FAST_Farm_IO.f90 b/glue-codes/fast-farm/src/FAST_Farm_IO.f90 index 32fce6d53c..d5f0bfab34 100644 --- a/glue-codes/fast-farm/src/FAST_Farm_IO.f90 +++ b/glue-codes/fast-farm/src/FAST_Farm_IO.f90 @@ -938,6 +938,12 @@ SUBROUTINE Farm_ReadPrimaryFile( InputFile, p, WD_InitInp, AWAE_InitInp, OutList CALL ReadVarWDefault( UnIn, InputFile, AWAE_InitInp%WrDisDT, "WrDisDT", "The time between vtk outputs [must be a multiple of the low resolution time step]", p%DT_low, ErrStat2, ErrMsg2, UnEc); if (Failed()) return + !---------------------- AXIS-ALIGNED PLANE SLICES (extent-controlled) [Feature 2] --- + call ReadPlaneSlicesBlock( UnIn, InputFile, AWAE_InitInp, p%DT_low, UnEc, ErrStat2, ErrMsg2 ); if (Failed()) return + + !---------------------- TERRAIN-FOLLOWING SAMPLING [Feature 3] --- + call ReadTerrainSlicesBlock( UnIn, InputFile, AWAE_InitInp, p%DT_low, UnEc, ErrStat2, ErrMsg2 ); if (Failed()) return + !---------------------- OUTPUT -------------------------------------------------- CALL ReadCom( UnIn, InputFile, 'Section Header: Output', ErrStat2, ErrMsg2, UnEc ); if (Failed()) return CALL ReadVar( UnIn, InputFile, p%SumPrint, "SumPrint", "Print summary data to .sum? (flag)", ErrStat2, ErrMsg2, UnEc); if (Failed()) return @@ -1212,5 +1218,586 @@ SUBROUTINE Farm_ValidateInput( p, WD_InitInp, AWAE_InitInp, ErrStat, ErrMsg ) END SUBROUTINE Farm_ValidateInput +!---------------------------------------------------------------------------------------------------------------------------------- +!> Read the required '--- AXIS-ALIGNED PLANE SLICES (extent-controlled) ---' +!! block for Feature 2 (axis-aligned planar sampling via user-specified +!! origin/plane/extents). The block layout is: +!! +!! --- AXIS-ALIGNED PLANE SLICES (extent-controlled) --- +!! NumPlaneSlices +!! WrPlaneDT +!! SliceName origin(m) plane extent1(m) extent2(m) +!! (-) (m,m,m) (-) (m) (m) +!! "T1_0D" (0 -300 0) YZ 600 400 +!! "hubXY" (-500 -300 0) XY 2000 800 +!! +!! This section is mandatory. Set NumPlaneSlices=0 to disable output. +subroutine ReadPlaneSlicesBlock( UnIn, InputFile, AWAE_InitInp, DT_low, UnEc, ErrStat, ErrMsg ) + integer(IntKi), intent(in ) :: UnIn + character(*), intent(in ) :: InputFile + type(AWAE_InputFileType), intent(inout) :: AWAE_InitInp + real(DbKi), intent(in ) :: DT_low + integer(IntKi), intent(in ) :: UnEc + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + + integer(IntKi) :: ErrStat2, k, ios, lineLen + character(ErrMsgLen) :: ErrMsg2 + character(1024) :: line + character(64) :: name + real(ReKi) :: origin(3), normal(3), extent1, extent2 + integer(IntKi), parameter :: maxSlices = 99 + character(*), parameter :: RoutineName = 'ReadPlaneSlicesBlock' + + ErrStat = ErrID_None + ErrMsg = '' + + ! Default: feature disabled + AWAE_InitInp%NumPlaneSlices = 0 + AWAE_InitInp%WrPlaneDT = DT_low + + ! Read the required section header + call ReadCom( UnIn, InputFile, 'Section Header: Axis-Aligned Plane Slices', ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + + ! We have the section header. Read the two scalars. + call ReadVar( UnIn, InputFile, AWAE_InitInp%NumPlaneSlices, "NumPlaneSlices", & + "Number of axis-aligned extent-controlled planar slices (-) [0 to 99]", & + ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + + if ( AWAE_InitInp%NumPlaneSlices < 0 .or. AWAE_InitInp%NumPlaneSlices > maxSlices ) then + call SetErrStat( ErrID_Fatal, 'NumPlaneSlices must be in the range [0, 99].', ErrStat, ErrMsg, RoutineName ) + return + end if + + call ReadVarWDefault( UnIn, InputFile, AWAE_InitInp%WrPlaneDT, "WrPlaneDT", & + "Feature 2 sampling period (s); DEFAULT falls back to DT_Low", & + DT_low, ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + + ! Two column-header lines (always present, even when NumPlaneSlices=0) + call ReadCom( UnIn, InputFile, 'Plane slices column names', ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + call ReadCom( UnIn, InputFile, 'Plane slices column units', ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + + if ( AWAE_InitInp%NumPlaneSlices == 0 ) return + + ! Allocate flat arrays + call AllocAry( AWAE_InitInp%PlaneSliceName, AWAE_InitInp%NumPlaneSlices, 'PlaneSliceName', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( AWAE_InitInp%PlaneSliceOrigin, 3, AWAE_InitInp%NumPlaneSlices, 'PlaneSliceOrigin', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( AWAE_InitInp%PlaneSliceNormal, 3, AWAE_InitInp%NumPlaneSlices, 'PlaneSliceNormal', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( AWAE_InitInp%PlaneSliceExtent1, AWAE_InitInp%NumPlaneSlices, 'PlaneSliceExtent1', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( AWAE_InitInp%PlaneSliceExtent2, AWAE_InitInp%NumPlaneSlices, 'PlaneSliceExtent2', ErrStat2, ErrMsg2 ); if (Failed()) return + + do k = 1, AWAE_InitInp%NumPlaneSlices + call ReadLine( UnIn, '', line, lineLen, ios ) + if (ios /= 0) then + call SetErrStat( ErrID_Fatal, 'Failed to read plane slice line '//trim(Num2LStr(k))//'.', ErrStat, ErrMsg, RoutineName ) + return + end if + if (UnEc > 0) write(UnEc,'(A)') trim(line) + + ! Strip parentheses and commas so parenthesized vectors "(700 700 95)" parse + ! as space-separated tokens. + call StripDelims( line ) + + call ParsePlaneSliceLine( line, name, origin, normal, extent1, extent2, ErrStat2, ErrMsg2 ) + if (ErrStat2 /= ErrID_None) then + call SetErrStat( ErrID_Fatal, 'Slice '//trim(Num2LStr(k))//': '//trim(ErrMsg2), ErrStat, ErrMsg, RoutineName ) + return + end if + + AWAE_InitInp%PlaneSliceName(k) = name + AWAE_InitInp%PlaneSliceOrigin(:,k) = origin + AWAE_InitInp%PlaneSliceNormal(:,k) = normal + AWAE_InitInp%PlaneSliceExtent1(k) = extent1 + AWAE_InitInp%PlaneSliceExtent2(k) = extent2 + end do + +contains + logical function Failed() + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + end function + + !> Replace '(', ')', and ',' with spaces so parenthesised vectors parse + !! as space-separated tokens. + subroutine StripDelims( s ) + character(*), intent(inout) :: s + integer :: i + do i = 1, len_trim(s) + if ( s(i:i) == '(' .or. s(i:i) == ')' .or. s(i:i) == ',' ) s(i:i) = ' ' + end do + end subroutine + + !> Parse a de-delimited slice row: NAME O_X O_Y O_Z PLANE E1 E2 + !! (plus trailing comment). Name may be quoted with double quotes. + !! PLANE must be one of XY, YZ, XZ (case-insensitive). + subroutine ParsePlaneSliceLine( raw, sliceName, o, n, e1, e2, err, msg ) + character(*), intent(in ) :: raw + character(*), intent( out) :: sliceName + real(ReKi), intent( out) :: o(3), n(3) + real(ReKi), intent( out) :: e1, e2 + integer(IntKi), intent( out) :: err + character(*), intent( out) :: msg + + character(len(raw)) :: work + character(64) :: tok(6) + character(8) :: planeTok + integer :: i, ipos, nTok, ts, te, iosLocal + + work = raw + msg = '' + err = ErrID_None + + ! Trim inline comments after '!' or '#' + i = scan(work, '!#') + if (i > 0) work(i:) = ' ' + + ! Extract the quoted-or-unquoted name and blank it out + call ExtractQuotedName( work, sliceName ) + + ! Tokenise remaining string by whitespace + nTok = 0 + ipos = 1 + do while (nTok < 6) + ! Skip blanks / tabs + do while (ipos <= len_trim(work)) + if (work(ipos:ipos) /= ' ' .and. work(ipos:ipos) /= char(9)) exit + ipos = ipos + 1 + end do + if (ipos > len_trim(work)) exit + ts = ipos + ! Advance to end of token + do while (ipos <= len_trim(work)) + if (work(ipos:ipos) == ' ' .or. work(ipos:ipos) == char(9)) exit + ipos = ipos + 1 + end do + te = ipos - 1 + nTok = nTok + 1 + tok(nTok) = work(ts:te) + end do + + if (nTok < 6) then + err = ErrID_Fatal + msg = 'Expected 6 columns after SliceName: Origin_X, Origin_Y, Origin_Z, Plane, Extent1, Extent2.' + return + end if + + ! Convert origin tokens to reals + read(tok(1), *, iostat=iosLocal) o(1) + if (iosLocal == 0) read(tok(2), *, iostat=iosLocal) o(2) + if (iosLocal == 0) read(tok(3), *, iostat=iosLocal) o(3) + if (iosLocal /= 0) then + err = ErrID_Fatal + msg = 'Failed to parse Origin (3) as numeric values.' + return + end if + + ! Convert extent tokens to reals + read(tok(5), *, iostat=iosLocal) e1 + if (iosLocal == 0) read(tok(6), *, iostat=iosLocal) e2 + if (iosLocal /= 0) then + err = ErrID_Fatal + msg = 'Failed to parse Extent1, Extent2 as numeric values.' + return + end if + + ! Convert plane keyword to axis-aligned normal vector + planeTok = tok(4) + call Conv2UC( planeTok ) + select case ( trim(planeTok) ) + case ('XY') + n = (/ 0.0_ReKi, 0.0_ReKi, 1.0_ReKi /) + case ('YZ') + n = (/ 1.0_ReKi, 0.0_ReKi, 0.0_ReKi /) + case ('XZ') + n = (/ 0.0_ReKi, 1.0_ReKi, 0.0_ReKi /) + case default + err = ErrID_Fatal + msg = 'Plane keyword "'//trim(planeTok)//'" is invalid. Must be XY, YZ, or XZ.' + end select + end subroutine + + !> Extract the first quoted-or-unquoted name from a line and blank it out + !! so subsequent list-directed reads see only numeric tokens. + subroutine ExtractQuotedName( s, nameOut ) + character(*), intent(inout) :: s + character(*), intent( out) :: nameOut + integer :: q1, q2, ns, ne + + nameOut = '' + + q1 = index(s, '"') + if (q1 > 0) then + q2 = index(s(q1+1:), '"') + if (q2 > 0) then + nameOut = adjustl(s(q1+1:q1+q2-1)) + s(q1:q1+q2) = repeat(' ', q2+1) + return + end if + end if + + ! Unquoted: first whitespace-delimited token + ns = verify(s, ' '//char(9)) + if (ns == 0) return + ne = ns + do + if (ne > len(s)) exit + if (s(ne:ne) == ' ' .or. s(ne:ne) == char(9)) exit + ne = ne + 1 + end do + nameOut = adjustl(s(ns:ne-1)) + s(ns:ne-1) = repeat(' ', ne-ns) + end subroutine + +end subroutine ReadPlaneSlicesBlock + +!---------------------------------------------------------------------------------------------------------------------------------- +!> Read the required '--- TERRAIN-FOLLOWING SAMPLING ---' block for +!! Feature 3 (terrain-following point-cloud sampling). The block layout is: +!! +!! --- TERRAIN-FOLLOWING SAMPLING --- +!! NumTerrainSlices [0 to 99] +!! WrTerrainDT +!! SliceName Offsets(m) OffsetNormal SourceType FileName +!! (-) (m,list) (-,-,-|dflt) (STL|Point) (quoted) +!! "terr" 100 200 300 default STL "terrain.stl" +!! "hubHt" 140.0 (0 0 1) Point "hub_points.txt" +!! +!! This section is mandatory. Set NumTerrainSlices=0 to disable output. +subroutine ReadTerrainSlicesBlock( UnIn, InputFile, AWAE_InitInp, DT_low, UnEc, ErrStat, ErrMsg ) + integer(IntKi), intent(in ) :: UnIn + character(*), intent(in ) :: InputFile + type(AWAE_InputFileType), intent(inout) :: AWAE_InitInp + real(DbKi), intent(in ) :: DT_low + integer(IntKi), intent(in ) :: UnEc + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + + integer(IntKi) :: ErrStat2, k, ios, lineLen, nOff, totOff + character(ErrMsgLen) :: ErrMsg2 + character(1024) :: line + character(64) :: name, srcTag + character(1024) :: fileName + real(ReKi) :: normal(3), offsets(64) + integer(IntKi), parameter :: maxSlices = 99 + integer(IntKi), parameter :: maxOffsets = 64 + integer(IntKi) :: perSliceNOff(maxSlices) + real(ReKi) :: perSliceOffsets(maxOffsets, maxSlices) + real(ReKi) :: perSliceNormal(3, maxSlices) + character(64) :: perSliceName(maxSlices) + character(1024) :: perSliceFile(maxSlices) + integer(IntKi) :: perSliceSrc (maxSlices) + character(*), parameter :: RoutineName = 'ReadTerrainSlicesBlock' + + ErrStat = ErrID_None + ErrMsg = '' + + ! Default: feature disabled + AWAE_InitInp%NumTerrainSlices = 0 + AWAE_InitInp%WrTerrainDT = DT_low + + ! Read the required section header + call ReadCom( UnIn, InputFile, 'Section Header: Terrain-Following Sampling', ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + + call ReadVar( UnIn, InputFile, AWAE_InitInp%NumTerrainSlices, "NumTerrainSlices", & + "Number of terrain-following point-cloud slices (-) [0 to 99]", & + ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + + if ( AWAE_InitInp%NumTerrainSlices < 0 .or. AWAE_InitInp%NumTerrainSlices > maxSlices ) then + call SetErrStat( ErrID_Fatal, 'NumTerrainSlices must be in the range [0, 99].', ErrStat, ErrMsg, RoutineName ) + return + end if + + call ReadVarWDefault( UnIn, InputFile, AWAE_InitInp%WrTerrainDT, "WrTerrainDT", & + "Feature 3 sampling period (s); DEFAULT falls back to DT_Low", & + DT_low, ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + + ! Two column-header lines (always present, even when NumTerrainSlices=0) + call ReadCom( UnIn, InputFile, 'Terrain slices column names', ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + call ReadCom( UnIn, InputFile, 'Terrain slices column units', ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + + if ( AWAE_InitInp%NumTerrainSlices == 0 ) return + + perSliceNOff = 0 + totOff = 0 + do k = 1, AWAE_InitInp%NumTerrainSlices + call ReadLine( UnIn, '', line, lineLen, ios ) + if (ios /= 0) then + call SetErrStat( ErrID_Fatal, 'Failed to read terrain slice line '//trim(Num2LStr(k))//'.', & + ErrStat, ErrMsg, RoutineName ); return + end if + if (UnEc > 0) write(UnEc,'(A)') trim(line) + + call ParseTerrainSliceLine( line, name, offsets, nOff, normal, srcTag, fileName, ErrStat2, ErrMsg2 ) + if (ErrStat2 /= ErrID_None) then + call SetErrStat( ErrID_Fatal, 'Slice '//trim(Num2LStr(k))//': '//trim(ErrMsg2), ErrStat, ErrMsg, RoutineName ) + return + end if + + perSliceName(k) = name + perSliceNOff(k) = nOff + perSliceOffsets(1:nOff, k) = offsets(1:nOff) + perSliceNormal(:,k) = normal + perSliceFile(k) = fileName + totOff = totOff + nOff + + ! Classify SourceType + call ToLower(srcTag) + select case (trim(srcTag)) + case ('stl') + perSliceSrc(k) = 1 + case ('point','points','pointcloud','point-cloud','point_cloud') + perSliceSrc(k) = 2 + case default + call SetErrStat( ErrID_Fatal, 'Slice '//trim(Num2LStr(k))// & + ': SourceType must be "STL" or "Point"; got "'//trim(srcTag)//'".', & + ErrStat, ErrMsg, RoutineName ) + return + end select + end do + + ! Allocate flat arrays now that we know sizes. + call AllocAry( AWAE_InitInp%TerrainSliceName, AWAE_InitInp%NumTerrainSlices, 'TerrainSliceName', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( AWAE_InitInp%TerrainSliceSourceType, AWAE_InitInp%NumTerrainSlices, 'TerrainSliceSourceType', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( AWAE_InitInp%TerrainSliceFileName, AWAE_InitInp%NumTerrainSlices, 'TerrainSliceFileName', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( AWAE_InitInp%TerrainSliceOffsetNormal, 3, AWAE_InitInp%NumTerrainSlices, 'TerrainSliceOffsetNormal', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( AWAE_InitInp%TerrainSliceOffsetsPacked, totOff, 'TerrainSliceOffsetsPacked', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( AWAE_InitInp%TerrainSliceOffsetIdx, AWAE_InitInp%NumTerrainSlices+1, 'TerrainSliceOffsetIdx', ErrStat2, ErrMsg2 ); if (Failed()) return + + AWAE_InitInp%TerrainSliceOffsetIdx(1) = 0 + do k = 1, AWAE_InitInp%NumTerrainSlices + AWAE_InitInp%TerrainSliceName(k) = perSliceName(k) + AWAE_InitInp%TerrainSliceSourceType(k) = perSliceSrc(k) + AWAE_InitInp%TerrainSliceFileName(k) = perSliceFile(k) + AWAE_InitInp%TerrainSliceOffsetNormal(:,k) = perSliceNormal(:,k) + AWAE_InitInp%TerrainSliceOffsetIdx(k+1) = AWAE_InitInp%TerrainSliceOffsetIdx(k) + perSliceNOff(k) + if (perSliceNOff(k) > 0) then + AWAE_InitInp%TerrainSliceOffsetsPacked( AWAE_InitInp%TerrainSliceOffsetIdx(k)+1 : AWAE_InitInp%TerrainSliceOffsetIdx(k+1) ) = & + perSliceOffsets(1:perSliceNOff(k), k) + end if + + ! Resolve relative filename against the primary input file's path + if ( PathIsRelative(AWAE_InitInp%TerrainSliceFileName(k)) ) then + AWAE_InitInp%TerrainSliceFileName(k) = trim(GetPathFromInput(InputFile))//trim(AWAE_InitInp%TerrainSliceFileName(k)) + end if + end do + +contains + logical function Failed() + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + end function + + function GetPathFromInput(fname) result(pathStr) + character(*), intent(in) :: fname + character(1024) :: pathStr + character(1024) :: dummyBase + call GetPath( fname, pathStr, dummyBase ) + end function + + subroutine ToLower(s) + character(*), intent(inout) :: s + integer :: i, c + do i = 1, len_trim(s) + c = iachar(s(i:i)) + if (c >= iachar('A') .and. c <= iachar('Z')) s(i:i) = achar(c + 32) + end do + end subroutine + + !> Parse a single terrain slice row after de-delimiting parens/commas. + !! Row format: NAME Offset1 [Offset2 ...] (Nx Ny Nz | default) SourceType "FileName" + !! The row is read once into a de-delimited string; we scan tokens with a + !! small state machine (tokens between the name and the (normal) group are + !! offsets; then 3 more numeric tokens are the normal or the 'default' + !! keyword; then SourceType keyword; then quoted filename). + subroutine ParseTerrainSliceLine( raw, sliceName, offList, nOffLocal, nrm, srcTagLocal, fnameLocal, err, msg ) + character(*), intent(in ) :: raw + character(*), intent( out) :: sliceName + real(ReKi), intent( out) :: offList(:) + integer(IntKi), intent( out) :: nOffLocal + real(ReKi), intent( out) :: nrm(3) + character(*), intent( out) :: srcTagLocal + character(*), intent( out) :: fnameLocal + integer(IntKi), intent( out) :: err + character(*), intent( out) :: msg + + character(len(raw)) :: work + character(1024) :: tok + integer :: tokCount, i, defPos, iosLocal + real(ReKi) :: rval + logical :: normalIsDefault + + err = ErrID_None + msg = '' + nrm = 0.0_ReKi + offList = 0.0_ReKi + nOffLocal = 0 + normalIsDefault = .false. + + work = raw + i = scan(work, '!#') + if (i > 0) work(i:) = ' ' + + ! SliceName must be quoted to distinguish it from the quoted FileName. + i = verify(work, ' '//char(9)) + if (i == 0 .or. work(i:i) /= '"') then + err = ErrID_Fatal; msg = 'SliceName must be a double-quoted string.' + return + end if + call ExtractQuoted( work, sliceName ) + if (len_trim(sliceName) == 0) then + err = ErrID_Fatal; msg = 'Missing slice name.' + return + end if + + call ExtractQuoted( work, fnameLocal ) + if (len_trim(fnameLocal) == 0) then + err = ErrID_Fatal; msg = 'Missing quoted FileName in terrain slice row.' + return + end if + + ! Strip parens/commas so vectors "(1 0 0)" parse as space-separated tokens + do i = 1, len_trim(work) + if ( work(i:i) == '(' .or. work(i:i) == ')' .or. work(i:i) == ',' ) work(i:i) = ' ' + end do + + ! Locate the case-insensitive "default" keyword if present. It replaces + ! the three normal tokens with an implicit signal that we should defer to + ! the STL per-facet normal. Blank out those chars so numeric scanning is + ! simpler afterwards. + defPos = FindKeywordCI( work, 'default' ) + if (defPos > 0) then + normalIsDefault = .true. + work(defPos:defPos+6) = ' ' + end if + + ! Now tokenise remaining words. Numeric tokens fill offList (variable + ! count). Non-numeric tokens are the SourceType keyword. If normal isn't + ! default, the last 3 numeric tokens before the SourceType keyword are + ! the normal vector. + tokCount = 0 + srcTagLocal = '' + do + call PopToken( work, tok ) + if (len_trim(tok) == 0) exit + + read(tok, *, iostat=iosLocal) rval + if (iosLocal == 0) then + ! Numeric: offset or normal component + if (tokCount >= size(offList)) then + err = ErrID_Fatal; msg = 'Too many offsets on a single line.' + return + end if + tokCount = tokCount + 1 + offList(tokCount) = rval + else + ! Non-numeric: SourceType keyword + srcTagLocal = adjustl(tok) + exit + end if + end do + + if (len_trim(srcTagLocal) == 0) then + err = ErrID_Fatal; msg = 'Missing SourceType (STL or Point).' + return + end if + + if (normalIsDefault) then + nOffLocal = tokCount + nrm = 0.0_ReKi ! sentinel: use per-facet normal from STL + else + if (tokCount < 3) then + err = ErrID_Fatal; msg = 'Need at least the (Nx Ny Nz) normal or the keyword "default".' + return + end if + nrm(1) = offList(tokCount-2) + nrm(2) = offList(tokCount-1) + nrm(3) = offList(tokCount ) + nOffLocal = tokCount - 3 + end if + + end subroutine + + subroutine ExtractQuoted( s, nameOut ) + character(*), intent(inout) :: s + character(*), intent( out) :: nameOut + integer :: q1, q2 + nameOut = '' + q1 = index(s, '"') + if (q1 == 0) return + q2 = index(s(q1+1:), '"') + if (q2 == 0) return + nameOut = adjustl(s(q1+1:q1+q2-1)) + s(q1:q1+q2) = repeat(' ', q2+1) + end subroutine + + subroutine ExtractBareToken( s, tokOut ) + character(*), intent(inout) :: s + character(*), intent( out) :: tokOut + integer :: ns, ne + tokOut = '' + ns = verify(s, ' '//char(9)) + if (ns == 0) return + ne = ns + do + if (ne > len(s)) exit + if (s(ne:ne) == ' ' .or. s(ne:ne) == char(9)) exit + ne = ne + 1 + end do + tokOut = adjustl(s(ns:ne-1)) + s(ns:ne-1) = repeat(' ', ne-ns) + end subroutine + + subroutine PopToken( s, tokOut ) + character(*), intent(inout) :: s + character(*), intent( out) :: tokOut + integer :: ns, ne + tokOut = '' + ns = verify(s, ' '//char(9)) + if (ns == 0) return + ne = ns + do + if (ne > len(s)) exit + if (s(ne:ne) == ' ' .or. s(ne:ne) == char(9)) exit + ne = ne + 1 + end do + tokOut = adjustl(s(ns:ne-1)) + s(ns:ne-1) = repeat(' ', ne-ns) + end subroutine + + integer function FindKeywordCI( s, kw ) + character(*), intent(in) :: s, kw + integer :: i, ls, lk + character(len(kw)) :: sub + FindKeywordCI = 0 + ls = len_trim(s); lk = len_trim(kw) + if (lk == 0 .or. ls < lk) return + do i = 1, ls - lk + 1 + sub = s(i:i+lk-1) + call ToLowerInner(sub) + if (sub == kw) then + ! word boundary check + if (i > 1) then + if (s(i-1:i-1) /= ' ' .and. s(i-1:i-1) /= char(9)) cycle + end if + if (i+lk <= ls) then + if (s(i+lk:i+lk) /= ' ' .and. s(i+lk:i+lk) /= char(9)) cycle + end if + FindKeywordCI = i + return + end if + end do + end function + + subroutine ToLowerInner(s) + character(*), intent(inout) :: s + integer :: i, c + do i = 1, len_trim(s) + c = iachar(s(i:i)) + if (c >= iachar('A') .and. c <= iachar('Z')) s(i:i) = achar(c + 32) + end do + end subroutine + +end subroutine ReadTerrainSlicesBlock + end module FAST_Farm_IO diff --git a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 index 5a2eabde5c..07e8fe988d 100644 --- a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 +++ b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 @@ -211,7 +211,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 @@ -620,6 +624,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 fdb66b91b6..289484c3ce 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 0f8360cc54..aa6aafccd2 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 @@ -59,6 +60,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 @@ -106,8 +113,871 @@ subroutine ExtractSlice( sliceType, s, s0, szs, sz1, sz2, ds, V, slice) end subroutine ExtractSlice !---------------------------------------------------------------------------------------------------------------------------------- -!> This subroutine +!> Bounded variant of ExtractSlice. Interpolates a single axis-aligned slice of the +!! 3-D low-resolution field V at coordinate `s` (metres, along the "thin" axis +!! implied by sliceType), then copies only the subregion +!! [i_lo:i_hi, j_lo:j_hi] into the output array. Any in-plane node whose +!! interpolated velocity comes from a grid location outside `[0,szs-1]` along +!! the thin axis is silently clamped (identical behaviour to ExtractSlice for +!! that axis). Nodes whose *in-plane* index falls outside the parent low-res +!! grid must be masked with NaN by the caller before calling this routine \u2014 +!! this routine trusts the (i_lo:i_hi, j_lo:j_hi) window and does not enforce +!! it against the low-res domain. !! +!! Used by the axis-aligned planar sampling feature (`NumPlaneSlices`). +subroutine ExtractSliceSub( sliceType, s, s0, szs, sz1, sz2, i_lo, i_hi, j_lo, j_hi, ds, V, slice) + + integer(IntKi), intent(in ) :: sliceType !< Type of slice: XYSlice, YZSlice, XZSlice + real(ReKi), intent(in ) :: s !< data value in meters of the interpolant along the thin axis + real(ReKi), intent(in ) :: s0 !< origin value in meters of the interpolant along the thin axis + integer(IntKi), intent(in ) :: szs !< number of grid points along the thin axis of V + integer(IntKi), intent(in ) :: sz1 !< 1st in-plane dimension of the parent low-res grid + integer(IntKi), intent(in ) :: sz2 !< 2nd in-plane dimension of the parent low-res grid + integer(IntKi), intent(in ) :: i_lo, i_hi !< 0-based inclusive index bounds along the parent's 1st in-plane axis; must satisfy 0 <= i_lo <= i_hi <= sz1-1 + integer(IntKi), intent(in ) :: j_lo, j_hi !< 0-based inclusive index bounds along the parent's 2nd in-plane axis; must satisfy 0 <= j_lo <= j_hi <= sz2-1 + real(ReKi), intent(in ) :: ds !< grid spacing along the thin axis + real(SiKi), intent(in ) :: V(:,0:,0:,0:) !< parent low-res field (3, 0:nX-1, 0:nY-1, 0:nZ-1) + real(SiKi), intent(inout) :: slice(:,0:,0:) !< output subregion; must be sized (3, 0:i_hi-i_lo, 0:j_hi-j_lo) + + integer(IntKi) :: s_grid0,s_grid1,i,j + real(SiKi) :: s_grid, sd + + ! Compute the interpolation cell along the thin axis (identical to ExtractSlice) + s_grid = real((s-s0)/ds,SiKi) + s_grid0 = floor(s_grid) + s_grid1 = s_grid0 + 1 + sd = (s_grid-real(s_grid0,SiKi)) + if (s_grid0 == (szs-1)) s_grid1 = s_grid0 + + ! Silently no-op if the bounds are inverted (allows callers to pass fully-out-of-domain slices) + if (i_hi < i_lo .or. j_hi < j_lo) return + + do j = j_lo, j_hi + do i = i_lo, i_hi + select case (sliceType) + case (XYSlice) + slice(:,i-i_lo,j-j_lo) = V(:,i,j,s_grid0)*(1.0_SiKi-sd) + V(:,i,j,s_grid1)*sd + case (YZSlice) + slice(:,i-i_lo,j-j_lo) = V(:,s_grid0,i,j)*(1.0_SiKi-sd) + V(:,s_grid1,i,j)*sd + case (XZSlice) + slice(:,i-i_lo,j-j_lo) = V(:,i,s_grid0,j)*(1.0_SiKi-sd) + V(:,i,s_grid1,j)*sd + end select + end do + end do + +end subroutine ExtractSliceSub +!---------------------------------------------------------------------------------------------------------------------------------- +!> Validate the Feature 2 (axis-aligned planar sampling) input records, derive +!! per-slice `NormalAxis` / `Axis1` / `Axis2` / `N1` / `N2` and the requested +!! parent-grid index window `(ILoReq, IHiReq, JLoReq, JHiReq)`. Off-axis normals +!! are fatal per plan §3.7. Out-of-domain windows are warnings only \u2014 the emit +!! step will fill any node whose parent index is outside `[0, n-1]` with NaN. +!! Also opens each valid slice's `.vts.series` sidecar for lifetime of the run. +subroutine ValidatePlaneSlices( p, InputFileData, errStat, errMsg ) + type(AWAE_ParameterType), intent(inout) :: p + type(AWAE_InputFileType), intent(in ) :: InputFileData + integer(IntKi), intent(inout) :: errStat + character(*), intent(inout) :: errMsg + + integer(IntKi) :: k, iaxis, n1, n2, iLo, iHi, jLo, jHi + real(ReKi) :: nrm(3), o(3), e1, e2, d1, d2 + real(ReKi), parameter :: nrmTol = 1.0e-6_ReKi + character(1024) :: seriesFile + integer(IntKi) :: seriesUn, ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + character(*), parameter :: RoutineName = 'ValidatePlaneSlices' + + if (p%NumPlaneSlices <= 0) return + + do k = 1, p%NumPlaneSlices + nrm = InputFileData%PlaneSliceNormal(:,k) + o = InputFileData%PlaneSliceOrigin(:,k) + e1 = InputFileData%PlaneSliceExtent1(k) + e2 = InputFileData%PlaneSliceExtent2(k) + + ! Axis-aligned normal check (fatal per plan §3.7) + if ( abs(nrm(1)-1.0_ReKi) < nrmTol .and. abs(nrm(2)) < nrmTol .and. abs(nrm(3)) < nrmTol ) then + iaxis = 1 + else if ( abs(nrm(2)-1.0_ReKi) < nrmTol .and. abs(nrm(1)) < nrmTol .and. abs(nrm(3)) < nrmTol ) then + iaxis = 2 + else if ( abs(nrm(3)-1.0_ReKi) < nrmTol .and. abs(nrm(1)) < nrmTol .and. abs(nrm(2)) < nrmTol ) then + iaxis = 3 + else + call SetErrStat( ErrID_Fatal, 'In slice "'//trim(InputFileData%PlaneSliceName(k))// & + '", the internal plane normal ('//trim(Num2LStr(nrm(1)))//', '// & + trim(Num2LStr(nrm(2)))//', '//trim(Num2LStr(nrm(3)))// & + ') is not axis-aligned. Valid Plane keywords are XY, YZ, XZ.', & + errStat, errMsg, RoutineName ) + cycle + end if + + if ( e1 <= 0.0_ReKi .or. e2 <= 0.0_ReKi ) then + call SetErrStat( ErrID_Fatal, 'In slice "'//trim(InputFileData%PlaneSliceName(k))// & + '", both extent1 and extent2 must be positive.', errStat, errMsg, RoutineName ) + cycle + end if + + ! Axis-1 = first non-normal axis in cyclic (X,Y,Z) order; Axis-2 = second. + select case (iaxis) + case (1); p%PlaneSliceAxis1(k) = 2; p%PlaneSliceAxis2(k) = 3 + case (2); p%PlaneSliceAxis1(k) = 1; p%PlaneSliceAxis2(k) = 3 + case (3); p%PlaneSliceAxis1(k) = 1; p%PlaneSliceAxis2(k) = 2 + end select + p%PlaneSliceNormalAxis(k) = iaxis + + d1 = p%LowRes%dXYZ( p%PlaneSliceAxis1(k) ) + d2 = p%LowRes%dXYZ( p%PlaneSliceAxis2(k) ) + n1 = nint( e1 / d1 ) + 1 + n2 = nint( e2 / d2 ) + 1 + p%PlaneSliceN1(k) = n1 + p%PlaneSliceN2(k) = n2 + + ! Requested in-plane index window (0-based) into the parent low-res grid + iLo = nint( ( o(p%PlaneSliceAxis1(k)) - p%LowRes%oXYZ(p%PlaneSliceAxis1(k)) ) / d1 ) + iHi = iLo + n1 - 1 + jLo = nint( ( o(p%PlaneSliceAxis2(k)) - p%LowRes%oXYZ(p%PlaneSliceAxis2(k)) ) / d2 ) + jHi = jLo + n2 - 1 + p%PlaneSliceILoReq(k) = iLo + p%PlaneSliceIHiReq(k) = iHi + p%PlaneSliceJLoReq(k) = jLo + p%PlaneSliceJHiReq(k) = jHi + + ! Check the plane's thin-axis location falls inside the domain + if ( o(iaxis) < p%LowRes%oXYZ(iaxis) .or. & + o(iaxis) > p%LowRes%oXYZ(iaxis) + real(p%LowRes%nXYZ(iaxis)-1, ReKi) * p%LowRes%dXYZ(iaxis) ) then + call SetErrStat( ErrID_Warn, 'In slice "'//trim(InputFileData%PlaneSliceName(k))// & + '", the plane location along the normal axis is outside the low-resolution domain. '// & + 'This slice will be skipped.', errStat, errMsg, RoutineName ) + p%PlaneSliceValid(k) = .false. + cycle + end if + + ! In-plane out-of-domain warning (still valid, just NaN-filled) + if ( iLo < 0 .or. iHi > p%LowRes%nXYZ(p%PlaneSliceAxis1(k))-1 .or. & + jLo < 0 .or. jHi > p%LowRes%nXYZ(p%PlaneSliceAxis2(k))-1 ) then + call SetErrStat( ErrID_Warn, 'In slice "'//trim(InputFileData%PlaneSliceName(k))// & + '", the requested plane extents straddle the low-resolution domain boundary. '// & + 'Nodes outside the domain will be written as NaN.', errStat, errMsg, RoutineName ) + end if + + ! Open the per-slice .vts.series sidecar + seriesFile = trim(p%OutFileFFvtkRoot)//".Plane."//trim(InputFileData%PlaneSliceName(k))//".vts.series" + call VTK_Series_Open( seriesFile, seriesUn, ErrStat2, ErrMsg2 ) + if (ErrStat2 < AbortErrLev) then + p%PlaneSliceSeriesUn(k) = seriesUn + else + call SetErrStat( ErrStat2, ErrMsg2, errStat, errMsg, RoutineName ) + p%PlaneSliceValid(k) = .false. + end if + end do + +end subroutine ValidatePlaneSlices +!---------------------------------------------------------------------------------------------------------------------------------- +!> Emit VTK output for every Feature 2 plane slice at the current low-res step. +!! Called from `AWAE_CalcOutput` whenever `mod(n, p%WrPlaneSkp) == 0`. +!! Out-of-domain nodes are pre-filled with IEEE quiet NaN so ParaView masks +!! them automatically. +subroutine EmitPlaneSlices( t, n, p, m, errStat, errMsg ) + use, intrinsic :: ieee_arithmetic, only : ieee_value, ieee_quiet_nan + real(DbKi), intent(in ) :: t + integer(IntKi), intent(in ) :: n + type(AWAE_ParameterType), intent(in ) :: p + type(AWAE_MiscVarType), intent(inout) :: m + integer(IntKi), intent( out) :: errStat + character(*), intent( out) :: errMsg + + integer(IntKi) :: k, iaxis, ax1, ax2, n1, n2 + integer(IntKi) :: iLoR, iHiR, jLoR, jHiR, iLoE, iHiE, jLoE, jHiE + integer(IntKi) :: i, j, sType, szThin, sz1, sz2 + real(ReKi) :: o(3), d1, d2, dThin, s0Thin, sVal, xyz(3) + real(SiKi), allocatable :: buf(:,:,:), sub(:,:,:) + real(ReKi), allocatable :: pts(:,:,:) + real(SiKi) :: nan_val + character(1024) :: fileName, seriesDir, seriesEntry + character(20) :: Tstr + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + character(*), parameter :: RoutineName = 'EmitPlaneSlices' + logical :: isFirst + + errStat = ErrID_None + errMsg = '' + nan_val = ieee_value( 0.0_SiKi, ieee_quiet_nan ) + + if (p%NumPlaneSlices <= 0) return + + ! Skip duplicate emit at the same low-res step. FARM_InitialCO calls + ! AWAE_CalcOutput twice at t=0 (once for validation and once for the actual + ! initial output); a second call would overwrite the identical .vts file + ! and add a duplicate entry to the .vts.series sidecar. + if ( m%LastPlaneSliceN == n ) return + m%LastPlaneSliceN = n + + write(Tstr, '(i' // trim(Num2LStr(p%VTK_tWidth)) //'.'// trim(Num2LStr(p%VTK_tWidth)) // ')') n/p%WrDisSkp1 + + do k = 1, p%NumPlaneSlices + if (.not. p%PlaneSliceValid(k)) cycle + + iaxis = p%PlaneSliceNormalAxis(k) + ax1 = p%PlaneSliceAxis1(k) + ax2 = p%PlaneSliceAxis2(k) + n1 = p%PlaneSliceN1(k) + n2 = p%PlaneSliceN2(k) + iLoR = p%PlaneSliceILoReq(k) + iHiR = p%PlaneSliceIHiReq(k) + jLoR = p%PlaneSliceJLoReq(k) + jHiR = p%PlaneSliceJHiReq(k) + o = p%PlaneSliceOrigin(:,k) + + d1 = p%LowRes%dXYZ(ax1) + d2 = p%LowRes%dXYZ(ax2) + dThin = p%LowRes%dXYZ(iaxis) + s0Thin = p%LowRes%oXYZ(iaxis) + sVal = o(iaxis) + szThin = p%LowRes%nXYZ(iaxis) + sz1 = p%LowRes%nXYZ(ax1) + sz2 = p%LowRes%nXYZ(ax2) + + select case (iaxis) + case (1); sType = YZSlice + case (2); sType = XZSlice + case (3); sType = XYSlice + end select + + ! Build the point cloud in the FAST.Farm global frame + allocate( pts(3, n1, n2), stat=ErrStat2 ) + if (ErrStat2 /= 0) then + call SetErrStat( ErrID_Fatal, 'Alloc pts failed', errStat, errMsg, RoutineName ); return + end if + do j = 1, n2 + do i = 1, n1 + xyz = o + xyz(ax1) = o(ax1) + real(i-1, ReKi) * d1 + xyz(ax2) = o(ax2) + real(j-1, ReKi) * d2 + pts(:,i,j) = xyz + end do + end do + + ! Allocate + NaN-fill the sample buffer + allocate( buf(3, 0:n1-1, 0:n2-1), stat=ErrStat2 ) + if (ErrStat2 /= 0) then + call SetErrStat( ErrID_Fatal, 'Alloc buf failed', errStat, errMsg, RoutineName ) + if (allocated(pts)) deallocate(pts) + return + end if + buf = nan_val + + ! Effective (in-domain) subregion + iLoE = max( 0, iLoR ) + iHiE = min( sz1-1, iHiR ) + jLoE = max( 0, jLoR ) + jHiE = min( sz2-1, jHiR ) + + if ( iHiE >= iLoE .and. jHiE >= jLoE ) then + allocate( sub(3, 0:iHiE-iLoE, 0:jHiE-jLoE), stat=ErrStat2 ) + if (ErrStat2 == 0) then + call ExtractSliceSub( sType, sVal, s0Thin, szThin, sz1, sz2, & + iLoE, iHiE, jLoE, jHiE, dThin, m%Vdist_low_full, sub ) + buf(:, iLoE-iLoR:iHiE-iLoR, jLoE-jLoR:jHiE-jLoR) = sub + deallocate(sub) + end if + end if + + ! Write .vts + append to .vts.series + fileName = trim(p%OutFileFFvtkRoot)//".Plane."//trim(p%PlaneSliceName(k))//"."//trim(Tstr)//".vts" + call WriteVTK_StructuredGrid_2D( fileName, & + "Plane slice "//trim(p%PlaneSliceName(k))//" at t="//trim(num2lstr(t))//"s", & + n1, n2, pts, reshape(buf, [3,n1,n2]), "Velocity", ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, errStat, errMsg, RoutineName ) + + if ( p%PlaneSliceSeriesUn(k) > 0 ) then + isFirst = ( m%PlaneSliceStepCount(k) == 0 ) + ! The series sidecar records file names relative to its own directory, + ! so strip the leading path component from fileName before appending. + call GetPath( fileName, seriesDir, seriesEntry ) + call VTK_Series_Append( p%PlaneSliceSeriesUn(k), trim(seriesEntry), t, isFirst, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, errStat, errMsg, RoutineName ) + m%PlaneSliceStepCount(k) = m%PlaneSliceStepCount(k) + 1 + end if + + deallocate(pts, buf) + end do + +end subroutine EmitPlaneSlices +!---------------------------------------------------------------------------------------------------------------------------------- +!> Read an ASCII or binary STL file and return the union of all vertex +!! coordinates plus a matching per-vertex facet normal (each vertex inherits +!! the normal of its owning facet; duplicate vertices are kept, so the point +!! count equals 3 * NumFacets). If the file cannot be opened the routine +!! returns ErrID_Fatal. The binary/ascii selection is inferred from the +!! first line: ASCII STL begins with the keyword "solid " followed by +!! printable text (binary STL headers may start with the byte pattern that +!! looks like "solid" too, so we cross-check by reading a few facet lines). +subroutine ReadSTLFile( fileName, Pts, Nrm, ErrStat, ErrMsg ) + character(*), intent(in ) :: fileName + real(ReKi), allocatable, intent( out) :: Pts(:,:) !< (3, 3*NumFacets) + real(ReKi), allocatable, intent( out) :: Nrm(:,:) !< (3, 3*NumFacets) - facet normal, per vertex + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + + integer(IntKi) :: Un, ios, ErrStat2, i, npts, nFacets, iface + character(ErrMsgLen) :: ErrMsg2 + character(1024) :: line + character(80) :: hdr + integer(kind=B4Ki) :: nFacets32 + real(SiKi) :: nrm32(3), v32(3) + integer(kind=B2Ki) :: attrByteCount + character(6) :: word6 + logical :: isAscii + character(*), parameter :: RoutineName = 'ReadSTLFile' + + ErrStat = ErrID_None + ErrMsg = '' + + !$OMP critical(fileopenNWTCio_critical) + call GetNewUnit(Un, ErrStat2, ErrMsg2) + ! Open as unformatted stream so we can inspect the header regardless of + ! whether the file is ASCII or binary. + open( unit=Un, file=trim(fileName), status='old', access='stream', & + form='unformatted', action='read', iostat=ios ) + !$OMP end critical(fileopenNWTCio_critical) + if (ios /= 0) then + call SetErrStat( ErrID_Fatal, 'Cannot open STL file "'//trim(fileName)//'"', ErrStat, ErrMsg, RoutineName ) + return + end if + + ! Read the first 6 characters. ASCII STLs start with "solid " and are + ! followed by a printable name + newline (\n or \r\n). Binary STLs also + ! typically start with "solid " in their 80-byte header but have no + ! newline in the first ~200 bytes. Test: read the first line; if it fits + ! neatly under 1024 chars and contains "facet " in the file's second + ! non-empty line, treat as ASCII. + read(Un, iostat=ios) word6 + close(Un) + isAscii = ( word6(1:5) == 'solid' .or. word6(1:5) == 'SOLID' ) + + if (isAscii) then + ! Reopen as text and probe for "facet" to distinguish from binary "solid" false positives. + !$OMP critical(fileopenNWTCio_critical) + open( unit=Un, file=trim(fileName), status='old', action='read', iostat=ios ) + !$OMP end critical(fileopenNWTCio_critical) + if (ios /= 0) then + call SetErrStat( ErrID_Fatal, 'Cannot reopen STL file "'//trim(fileName)//'"', ErrStat, ErrMsg, RoutineName ) + return + end if + + ! Count facets first + nFacets = 0 + do + read(Un, '(A)', iostat=ios) line + if (ios /= 0) exit + line = adjustl(line) + call SToLower(line) + if (index(line, 'facet normal') == 1) nFacets = nFacets + 1 + end do + + if (nFacets == 0) then + ! No "facet normal" found: the file was misidentified as ASCII or is empty. + !$OMP critical(fileopenNWTCio_critical) + close(Un) + !$OMP end critical(fileopenNWTCio_critical) + call ReadSTL_Binary( fileName, Pts, Nrm, ErrStat, ErrMsg ) + return + end if + + npts = 3 * nFacets + allocate(Pts(3, npts), Nrm(3, npts), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat( ErrID_Fatal, 'Alloc Pts/Nrm for '//trim(Num2LStr(npts))//' STL vertices failed', & + ErrStat, ErrMsg, RoutineName ) + return + end if + + rewind(Un) + iface = 0 + npts = 0 + do + read(Un, '(A)', iostat=ios) line + if (ios /= 0) exit + line = adjustl(line) + call SToLower(line) + if (index(line, 'facet normal') == 1) then + iface = iface + 1 + read(line(len('facet normal')+1:), *, iostat=ios) nrm32(1), nrm32(2), nrm32(3) + if (ios /= 0) then + call SetErrStat( ErrID_Fatal, 'Failed to parse normal for facet '//trim(Num2LStr(iface))//' in "'//trim(fileName)//'".', ErrStat, ErrMsg, RoutineName ) + close(Un); return + end if + else if (index(line, 'vertex') == 1) then + read(line(len('vertex')+1:), *, iostat=ios) v32(1), v32(2), v32(3) + if (ios /= 0) then + call SetErrStat( ErrID_Fatal, 'Failed to parse vertex in "'//trim(fileName)//'".', ErrStat, ErrMsg, RoutineName ) + close(Un); return + end if + npts = npts + 1 + Pts(:, npts) = real(v32, ReKi) + Nrm(:, npts) = real(nrm32, ReKi) + end if + end do + close(Un) + + else + call ReadSTL_Binary( fileName, Pts, Nrm, ErrStat, ErrMsg ) + end if + +contains + subroutine SToLower(s) + character(*), intent(inout) :: s + integer :: kk, cc + do kk = 1, len_trim(s) + cc = iachar(s(kk:kk)) + if (cc >= iachar('A') .and. cc <= iachar('Z')) s(kk:kk) = achar(cc + 32) + end do + end subroutine + + subroutine ReadSTL_Binary( fname, Pts_, Nrm_, es, em ) + character(*), intent(in ) :: fname + real(ReKi), allocatable, intent(out) :: Pts_(:,:), Nrm_(:,:) + integer(IntKi), intent(out) :: es + character(*), intent(out) :: em + integer :: Un2, iosB, iB, jB, kB + integer(kind=B4Ki) :: nFB + real(SiKi) :: nrmB(3), vB(3) + integer(kind=B2Ki) :: attrB + + es = ErrID_None; em = '' + + !$OMP critical(fileopenNWTCio_critical) + call GetNewUnit(Un2, es, em) + open( unit=Un2, file=trim(fname), status='old', access='stream', & + form='unformatted', action='read', iostat=iosB ) + !$OMP end critical(fileopenNWTCio_critical) + if (iosB /= 0) then + es = ErrID_Fatal + em = 'Cannot open binary STL "'//trim(fname)//'"' + return + end if + + read(Un2, iostat=iosB) hdr + read(Un2, iostat=iosB) nFB + if (iosB /= 0 .or. nFB <= 0) then + close(Un2) + es = ErrID_Fatal + em = 'Failed to read binary STL header/facet count from "'//trim(fname)//'"' + return + end if + + allocate(Pts_(3, 3*int(nFB,IntKi)), Nrm_(3, 3*int(nFB,IntKi)), stat=iosB) + if (iosB /= 0) then + close(Un2) + es = ErrID_Fatal + em = 'Alloc for binary STL vertices failed' + return + end if + + kB = 0 + do iB = 1, int(nFB, IntKi) + read(Un2, iostat=iosB) nrmB + if (iosB /= 0) exit + do jB = 1, 3 + read(Un2, iostat=iosB) vB + if (iosB /= 0) exit + kB = kB + 1 + Pts_(:, kB) = real(vB, ReKi) + Nrm_(:, kB) = real(nrmB, ReKi) + end do + read(Un2, iostat=iosB) attrB + end do + close(Un2) + end subroutine + +end subroutine ReadSTLFile +!---------------------------------------------------------------------------------------------------------------------------------- +!> Read a plain-text or CSV point-cloud file: one point per line, three +!! numeric columns `x y z` in the FAST.Farm global coordinate frame. White +!! space OR comma delimiters; leading '!' or '#' lines are comments; an +!! optional header row is accepted as long as it fails numeric parsing and +!! appears before the first data row. Point count is derived from the file. +subroutine ReadPointCloudFile( fileName, Pts, ErrStat, ErrMsg ) + character(*), intent(in ) :: fileName + real(ReKi), allocatable, intent( out) :: Pts(:,:) + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + + integer(IntKi) :: Un, ios, ErrStat2, npts, k + character(ErrMsgLen) :: ErrMsg2 + character(1024) :: line + real(ReKi) :: x, y, z + logical :: sawFirstData + character(*), parameter :: RoutineName = 'ReadPointCloudFile' + + ErrStat = ErrID_None + ErrMsg = '' + + !$OMP critical(fileopenNWTCio_critical) + call GetNewUnit(Un, ErrStat2, ErrMsg2) + open( unit=Un, file=trim(fileName), status='old', action='read', iostat=ios ) + !$OMP end critical(fileopenNWTCio_critical) + if (ios /= 0) then + call SetErrStat( ErrID_Fatal, 'Cannot open point-cloud file "'//trim(fileName)//'"', ErrStat, ErrMsg, RoutineName ) + return + end if + + ! First pass: count data lines + npts = 0 + sawFirstData = .false. + do + read(Un, '(A)', iostat=ios) line + if (ios /= 0) exit + call NormalizeLine( line ) + if (len_trim(line) == 0) cycle + read(line, *, iostat=ios) x, y, z + if (ios == 0) then + npts = npts + 1 + sawFirstData = .true. + else + if (sawFirstData) then + call SetErrStat( ErrID_Fatal, 'Malformed row in point cloud "'//trim(fileName)//'": '//trim(line), & + ErrStat, ErrMsg, RoutineName ) + close(Un); return + end if + ! Header row before any data: silently skip + end if + end do + + if (npts == 0) then + close(Un) + call SetErrStat( ErrID_Warn, 'Point-cloud file "'//trim(fileName)//'" contained zero data rows.', & + ErrStat, ErrMsg, RoutineName ) + return + end if + + allocate(Pts(3, npts), stat=ErrStat2) + if (ErrStat2 /= 0) then + close(Un) + call SetErrStat( ErrID_Fatal, 'Alloc Pts for '//trim(Num2LStr(npts))//' point-cloud vertices failed', & + ErrStat, ErrMsg, RoutineName ) + return + end if + + ! Second pass: fill + rewind(Un) + k = 0 + do + read(Un, '(A)', iostat=ios) line + if (ios /= 0) exit + call NormalizeLine( line ) + if (len_trim(line) == 0) cycle + read(line, *, iostat=ios) x, y, z + if (ios /= 0) cycle + k = k + 1 + Pts(1,k) = x + Pts(2,k) = y + Pts(3,k) = z + end do + close(Un) + +contains + !> Strip inline comments, replace commas with spaces so both plain-text + !! and CSV files parse under list-directed reads. + subroutine NormalizeLine( s ) + character(*), intent(inout) :: s + integer :: kk + kk = scan(s, '!#') + if (kk > 0) s(kk:) = ' ' + do kk = 1, len_trim(s) + if (s(kk:kk) == ',') s(kk:kk) = ' ' + end do + end subroutine + +end subroutine ReadPointCloudFile +!---------------------------------------------------------------------------------------------------------------------------------- +!> Read every terrain slice's source file, expand each source with its list of +!! offsets, and pack the resulting point coordinates into `p%TerrainSlicePtsFlat`. +!! Emits an `ErrID_Info` size advisory (per plan §2.2.2) when a slice exceeds +!! ~10 M sample locations. Opens each slice's `.vtp.series` sidecar for the +!! lifetime of the run. Also emits an `ErrID_Warn` at init when part of the +!! slice sits outside the low-res domain (out-of-domain points are NaN-filled +!! at emit time). +subroutine ValidateTerrainSlices( p, InputFileData, ErrStat, ErrMsg ) + type(AWAE_ParameterType), intent(inout) :: p + type(AWAE_InputFileType), intent(in ) :: InputFileData + integer(IntKi), intent(inout) :: ErrStat + character(*), intent(inout) :: ErrMsg + + integer(IntKi) :: k, oStart, oEnd, nOff, iOff, iSrc, iOut, totOut + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + real(ReKi), allocatable :: srcPts(:,:), srcNrm(:,:) + real(ReKi) :: gNrm(3), o + character(1024) :: seriesFile + integer(IntKi) :: seriesUn + integer(IntKi) :: nSrc, nOutside + real(ReKi) :: xyz(3) + character(*), parameter :: RoutineName = 'ValidateTerrainSlices' + integer(IntKi), parameter :: sizeAdvisoryThreshold = 10000000 + + if (p%NumTerrainSlices <= 0) return + + ! Pass 1: read each source file, compute total point count per slice + call AllocAry( p%TerrainSliceNPtsTotal, p%NumTerrainSlices, 'p%TerrainSliceNPtsTotal', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%TerrainSlicePtsIdx, p%NumTerrainSlices+1, 'p%TerrainSlicePtsIdx', ErrStat2, ErrMsg2 ); if (Failed()) return + + p%TerrainSlicePtsIdx(1) = 0 + totOut = 0 + do k = 1, p%NumTerrainSlices + oStart = InputFileData%TerrainSliceOffsetIdx(k) + oEnd = InputFileData%TerrainSliceOffsetIdx(k+1) + nOff = oEnd - oStart + if (nOff <= 0) then + call SetErrStat( ErrID_Warn, 'Terrain slice "'//trim(InputFileData%TerrainSliceName(k))// & + '" has zero offsets; a single sheet at offset=0 will be produced.', ErrStat, ErrMsg, RoutineName ) + nOff = 1 ! implicit offset 0 + end if + + ! Read source file to count vertices + if (InputFileData%TerrainSliceSourceType(k) == 1) then + call ReadSTLFile( InputFileData%TerrainSliceFileName(k), srcPts, srcNrm, ErrStat2, ErrMsg2 ) + else + call ReadPointCloudFile( InputFileData%TerrainSliceFileName(k), srcPts, ErrStat2, ErrMsg2 ) + end if + if (Failed()) return + + nSrc = 0 + if (allocated(srcPts)) nSrc = size(srcPts, 2) + + p%TerrainSliceNPtsTotal(k) = nSrc * nOff + p%TerrainSlicePtsIdx(k+1) = p%TerrainSlicePtsIdx(k) + p%TerrainSliceNPtsTotal(k) + totOut = totOut + p%TerrainSliceNPtsTotal(k) + + if ( p%TerrainSliceNPtsTotal(k) >= sizeAdvisoryThreshold ) then + call SetErrStat( ErrID_Info, 'Terrain slice "'//trim(InputFileData%TerrainSliceName(k))// & + '" has '//trim(Num2LStr(p%TerrainSliceNPtsTotal(k)))// & + ' sample locations (~'//trim(Num2LStr(real(p%TerrainSliceNPtsTotal(k)*28)/1e6))// & + ' MB persistent RAM).', ErrStat, ErrMsg, RoutineName ) + end if + + if (allocated(srcPts)) deallocate(srcPts) + if (allocated(srcNrm)) deallocate(srcNrm) + end do + + ! Pass 2: allocate the flat coordinate array + call AllocAry( p%TerrainSlicePtsFlat, 3, max(totOut,1_IntKi), 'p%TerrainSlicePtsFlat', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%TerrainSliceName, p%NumTerrainSlices, 'p%TerrainSliceName', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%TerrainSliceValid, p%NumTerrainSlices, 'p%TerrainSliceValid', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%TerrainSliceSeriesUn, p%NumTerrainSlices, 'p%TerrainSliceSeriesUn', ErrStat2, ErrMsg2 ); if (Failed()) return + p%TerrainSliceValid = .true. + p%TerrainSliceSeriesUn = -1 + + ! Pass 3: re-read each file, expand with offsets, write into the flat array + do k = 1, p%NumTerrainSlices + p%TerrainSliceName(k) = InputFileData%TerrainSliceName(k) + + oStart = InputFileData%TerrainSliceOffsetIdx(k) + oEnd = InputFileData%TerrainSliceOffsetIdx(k+1) + nOff = max(1, oEnd - oStart) + + if (InputFileData%TerrainSliceSourceType(k) == 1) then + call ReadSTLFile( InputFileData%TerrainSliceFileName(k), srcPts, srcNrm, ErrStat2, ErrMsg2 ) + else + call ReadPointCloudFile( InputFileData%TerrainSliceFileName(k), srcPts, ErrStat2, ErrMsg2 ) + end if + if (Failed()) return + if (.not. allocated(srcPts)) cycle + nSrc = size(srcPts, 2) + if (nSrc == 0) cycle + + gNrm = InputFileData%TerrainSliceOffsetNormal(:,k) + + iOut = p%TerrainSlicePtsIdx(k) + nOutside = 0 + do iOff = 1, nOff + if (nOff == 1 .and. oEnd == oStart) then + o = 0.0_ReKi + else + o = InputFileData%TerrainSliceOffsetsPacked( oStart + iOff ) + end if + + do iSrc = 1, nSrc + iOut = iOut + 1 + if (all(abs(gNrm) < 1.0e-6_ReKi) .and. allocated(srcNrm)) then + ! Use per-facet normal from STL + xyz = srcPts(:,iSrc) + o * srcNrm(:,iSrc) + else + xyz = srcPts(:,iSrc) + o * gNrm + end if + p%TerrainSlicePtsFlat(:, iOut) = xyz + + if ( xyz(1) < p%LowRes%oXYZ(1) .or. xyz(1) > p%LowRes%oXYZ(1) + real(p%LowRes%nXYZ(1)-1,ReKi)*p%LowRes%dXYZ(1) .or. & + xyz(2) < p%LowRes%oXYZ(2) .or. xyz(2) > p%LowRes%oXYZ(2) + real(p%LowRes%nXYZ(2)-1,ReKi)*p%LowRes%dXYZ(2) .or. & + xyz(3) < p%LowRes%oXYZ(3) .or. xyz(3) > p%LowRes%oXYZ(3) + real(p%LowRes%nXYZ(3)-1,ReKi)*p%LowRes%dXYZ(3) ) then + nOutside = nOutside + 1 + end if + end do + end do + + if (nOutside > 0) then + call SetErrStat( ErrID_Warn, 'Terrain slice "'//trim(p%TerrainSliceName(k))//'": '// & + trim(Num2LStr(nOutside))//' of '//trim(Num2LStr(p%TerrainSliceNPtsTotal(k)))// & + ' sample points fall outside the low-resolution domain and will be written as NaN.', & + ErrStat, ErrMsg, RoutineName ) + end if + if (nOutside == p%TerrainSliceNPtsTotal(k)) then + call SetErrStat( ErrID_Warn, 'Terrain slice "'//trim(p%TerrainSliceName(k))// & + '" is entirely outside the low-resolution domain; skipping.', ErrStat, ErrMsg, RoutineName ) + p%TerrainSliceValid(k) = .false. + end if + + ! Open .vtp.series sidecar + if ( p%TerrainSliceValid(k) ) then + seriesFile = trim(p%OutFileFFvtkRoot)//".TerrSlice."//trim(p%TerrainSliceName(k))//".vtp.series" + call VTK_Series_Open( seriesFile, seriesUn, ErrStat2, ErrMsg2 ) + if (ErrStat2 < AbortErrLev) then + p%TerrainSliceSeriesUn(k) = seriesUn + else + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + p%TerrainSliceValid(k) = .false. + end if + end if + + if (allocated(srcPts)) deallocate(srcPts) + if (allocated(srcNrm)) deallocate(srcNrm) + end do + +contains + logical function Failed() + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + end function + +end subroutine ValidateTerrainSlices +!---------------------------------------------------------------------------------------------------------------------------------- +!> Emit VTK PolyData output for every Feature 3 terrain slice at the current +!! low-res step. Trilinear-samples m%Vdist_low_full at each pre-computed point; +!! out-of-domain points are set to IEEE quiet NaN (ParaView masks them). +subroutine EmitTerrainSlices( t, n, p, m, ErrStat, ErrMsg ) + use, intrinsic :: ieee_arithmetic, only : ieee_value, ieee_quiet_nan + real(DbKi), intent(in ) :: t + integer(IntKi), intent(in ) :: n + type(AWAE_ParameterType), intent(in ) :: p + type(AWAE_MiscVarType), intent(inout) :: m + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + + integer(IntKi) :: k, i, iStart, iEnd, nPts, ErrStat2 + real(SiKi), allocatable :: vel(:,:) + real(ReKi), allocatable :: ptsK(:,:) + real(SiKi) :: nan_val, u, v, w + real(ReKi) :: xr, yr, zr, fx, fy, fz + integer(IntKi) :: ix, iy, iz + character(ErrMsgLen) :: ErrMsg2 + character(1024) :: fileName, seriesDir, seriesEntry + character(20) :: Tstr + logical :: isFirst + character(*), parameter :: RoutineName = 'EmitTerrainSlices' + + ErrStat = ErrID_None + ErrMsg = '' + + if (p%NumTerrainSlices <= 0) return + + ! Guard the duplicate-emit at t=0 (same as Feature 2) + if ( m%LastTerrainSliceN == n ) return + m%LastTerrainSliceN = n + + nan_val = ieee_value( 0.0_SiKi, ieee_quiet_nan ) + + write(Tstr, '(i' // trim(Num2LStr(p%VTK_tWidth)) //'.'// trim(Num2LStr(p%VTK_tWidth)) // ')') n/p%WrDisSkp1 + + do k = 1, p%NumTerrainSlices + if (.not. p%TerrainSliceValid(k)) cycle + + iStart = p%TerrainSlicePtsIdx(k) + 1 + iEnd = p%TerrainSlicePtsIdx(k+1) + nPts = iEnd - iStart + 1 + if (nPts <= 0) cycle + + allocate(ptsK(3, nPts), vel(3, nPts), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat( ErrID_Fatal, 'Alloc temp Pts/Vel failed', ErrStat, ErrMsg, RoutineName ) + return + end if + ptsK = p%TerrainSlicePtsFlat(:, iStart:iEnd) + + do i = 1, nPts + ! Convert world coord to floating grid index + xr = (ptsK(1,i) - p%LowRes%oXYZ(1)) / p%LowRes%dXYZ(1) + yr = (ptsK(2,i) - p%LowRes%oXYZ(2)) / p%LowRes%dXYZ(2) + zr = (ptsK(3,i) - p%LowRes%oXYZ(3)) / p%LowRes%dXYZ(3) + + if ( xr < 0.0_ReKi .or. xr > real(p%LowRes%nXYZ(1)-1,ReKi) .or. & + yr < 0.0_ReKi .or. yr > real(p%LowRes%nXYZ(2)-1,ReKi) .or. & + zr < 0.0_ReKi .or. zr > real(p%LowRes%nXYZ(3)-1,ReKi) ) then + vel(:,i) = nan_val + cycle + end if + + ix = floor(xr); iy = floor(yr); iz = floor(zr) + if (ix > p%LowRes%nXYZ(1)-2) ix = p%LowRes%nXYZ(1)-2 + if (iy > p%LowRes%nXYZ(2)-2) iy = p%LowRes%nXYZ(2)-2 + if (iz > p%LowRes%nXYZ(3)-2) iz = p%LowRes%nXYZ(3)-2 + if (ix < 0) ix = 0 + if (iy < 0) iy = 0 + if (iz < 0) iz = 0 + fx = xr - real(ix,ReKi); fy = yr - real(iy,ReKi); fz = zr - real(iz,ReKi) + + u = trilerp_c( 1, ix, iy, iz, fx, fy, fz ) + v = trilerp_c( 2, ix, iy, iz, fx, fy, fz ) + w = trilerp_c( 3, ix, iy, iz, fx, fy, fz ) + vel(1,i) = u; vel(2,i) = v; vel(3,i) = w + end do + + fileName = trim(p%OutFileFFvtkRoot)//".TerrSlice."//trim(p%TerrainSliceName(k))//"."//trim(Tstr)//".vtp" + call WriteVTK_PolyData( fileName, & + "Terrain slice "//trim(p%TerrainSliceName(k))//" at t="//trim(num2lstr(t))//"s", & + ptsK, vel, "Velocity", ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + + if ( p%TerrainSliceSeriesUn(k) > 0 ) then + isFirst = ( m%TerrainSliceStepCount(k) == 0 ) + call GetPath( fileName, seriesDir, seriesEntry ) + call VTK_Series_Append( p%TerrainSliceSeriesUn(k), trim(seriesEntry), t, isFirst, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + m%TerrainSliceStepCount(k) = m%TerrainSliceStepCount(k) + 1 + end if + + deallocate(ptsK, vel) + end do + +contains + real(SiKi) function trilerp_c( icomp, ix0, iy0, iz0, fxL, fyL, fzL ) result(val) + integer(IntKi), intent(in) :: icomp, ix0, iy0, iz0 + real(ReKi), intent(in) :: fxL, fyL, fzL + real(SiKi) :: c00, c01, c10, c11, c0, c1 + c00 = m%Vdist_low_full(icomp, ix0, iy0, iz0) * (1.0_ReKi - fxL) + & + m%Vdist_low_full(icomp, ix0+1, iy0, iz0) * fxL + c01 = m%Vdist_low_full(icomp, ix0, iy0, iz0+1) * (1.0_ReKi - fxL) + & + m%Vdist_low_full(icomp, ix0+1, iy0, iz0+1) * fxL + c10 = m%Vdist_low_full(icomp, ix0, iy0+1, iz0) * (1.0_ReKi - fxL) + & + m%Vdist_low_full(icomp, ix0+1, iy0+1, iz0) * fxL + c11 = m%Vdist_low_full(icomp, ix0, iy0+1, iz0+1) * (1.0_ReKi - fxL) + & + m%Vdist_low_full(icomp, ix0+1, iy0+1, iz0+1) * fxL + c0 = c00 * (1.0_ReKi - fyL) + c10 * fyL + c1 = c01 * (1.0_ReKi - fyL) + c11 * fyL + val = c0 * (1.0_ReKi - fzL) + c1 * fzL + end function +end subroutine EmitTerrainSlices +!---------------------------------------------------------------------------------------------------------------------------------- +!> 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 @@ -816,8 +1686,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 @@ -877,97 +1747,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 @@ -1172,7 +2017,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 @@ -1251,15 +2096,69 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO p%OutDisWindX = InitInp%InputFileData%OutDisWindX p%OutDisWindY = InitInp%InputFileData%OutDisWindY + ! --- Feature 2: axis-aligned planar sampling (extent-controlled) + p%NumPlaneSlices = InitInp%InputFileData%NumPlaneSlices + p%WrPlaneSkp = max( 1_IntKi, nint( InitInp%InputFileData%WrPlaneDT / p%dt_low ) ) + if (p%NumPlaneSlices > 0) then + call AllocAry( p%PlaneSliceName, p%NumPlaneSlices, 'p%PlaneSliceName', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceOrigin, 3, p%NumPlaneSlices, 'p%PlaneSliceOrigin', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceNormalAxis, p%NumPlaneSlices, 'p%PlaneSliceNormalAxis', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceAxis1, p%NumPlaneSlices, 'p%PlaneSliceAxis1', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceAxis2, p%NumPlaneSlices, 'p%PlaneSliceAxis2', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceN1, p%NumPlaneSlices, 'p%PlaneSliceN1', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceN2, p%NumPlaneSlices, 'p%PlaneSliceN2', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceILoReq, p%NumPlaneSlices, 'p%PlaneSliceILoReq', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceIHiReq, p%NumPlaneSlices, 'p%PlaneSliceIHiReq', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceJLoReq, p%NumPlaneSlices, 'p%PlaneSliceJLoReq', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceJHiReq, p%NumPlaneSlices, 'p%PlaneSliceJHiReq', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceValid, p%NumPlaneSlices, 'p%PlaneSliceValid', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%PlaneSliceSeriesUn, p%NumPlaneSlices, 'p%PlaneSliceSeriesUn', ErrStat2, ErrMsg2 ); if (Failed()) return + p%PlaneSliceName = InitInp%InputFileData%PlaneSliceName + p%PlaneSliceOrigin = InitInp%InputFileData%PlaneSliceOrigin + p%PlaneSliceValid = .true. + p%PlaneSliceSeriesUn = -1 + end if + + ! --- Feature 3: terrain-following point-cloud sampling + p%NumTerrainSlices = InitInp%InputFileData%NumTerrainSlices + p%WrTerrainSkp = max( 1_IntKi, nint( InitInp%InputFileData%WrTerrainDT / p%dt_low ) ) + ! The rest of the terrain-slice parameter arrays are allocated inside + ! ValidateTerrainSlices (see below), because their sizes depend on reading + ! the per-slice source files. + ! --- 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) + ! placeholder for writing planes -- this will eventually be an input (revise logic here then) + p%WrPlanes = .true. + else if (p%NumPlaneSlices>0 .or. p%NumTerrainSlices>0) then + ! Features 2/3 (NumPlaneSlices/NumTerrainSlices) are independent of the wake-plane / disturbed-wind + ! output; still need the vtk_ff directory but must NOT enable WrPlanes. + call MKDIR(OutFileVTKDir) end if + ! 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; @@ -1414,6 +2313,23 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO end do if (errStat >= AbortErrLev) return + ! Feature 2: axis-aligned planar sampling + ! Validate each slice's normal (fatal on off-axis), derive its NormalAxis/Axis1/Axis2 + ! and (N1, N2) grid dimensions from the low-res spacing, and pre-compute the + ! requested (i_lo, i_hi, j_lo, j_hi) index window into the parent low-res grid. + ! Warn if the requested window straddles the domain boundary (out-of-domain + ! nodes will be written as NaN at emit time). + call ValidatePlaneSlices( p, InitInp%InputFileData, errStat, errMsg ) + if (errStat >= AbortErrLev) return + + ! Feature 3: terrain-following point-cloud sampling + ! Read each STL/point-cloud source file, expand with the per-slice offsets, + ! and pack the resulting point coordinates into p%TerrainSlicePtsFlat. + ! Opens each slice's .vtp.series sidecar and warns about out-of-domain + ! points (NaN-filled at emit time). + call ValidateTerrainSlices( p, InitInp%InputFileData, errStat, errMsg ) + if (errStat >= AbortErrLev) return + !---------------------------------------------------------------------------- ! Initialize inputs @@ -1428,10 +2344,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 !---------------------------------------------------------------------------- @@ -1457,9 +2377,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 @@ -1471,33 +2391,51 @@ 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 + if ( p%NumPlaneSlices > 0 ) then + call AllocAry( m%PlaneSliceStepCount, p%NumPlaneSlices, 'm%PlaneSliceStepCount', ErrStat2, ErrMsg2 ); if (Failed()) return + m%PlaneSliceStepCount = 0 + end if + m%LastPlaneSliceN = -1 + if ( p%NumTerrainSlices > 0 ) then + call AllocAry( m%TerrainSliceStepCount, p%NumTerrainSlices, 'm%TerrainSliceStepCount', ErrStat2, ErrMsg2 ); if (Failed()) return + m%TerrainSliceStepCount = 0 + end if + m%LastTerrainSliceN = -1 ! miscvars to avoid the allocation per timestep allocate(m%Vamb_lowpol( 3, 0:p%n_rp_max*8 ), STAT=errStat2); if (Failed0('m%Vamb_lowpol.' )) return; 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; @@ -1526,9 +2464,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 @@ -1660,6 +2598,44 @@ 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 + + ! Feature 2: close any open .vts.series sidecars + if (allocated(p%PlaneSliceSeriesUn)) then + do nt = 1, size(p%PlaneSliceSeriesUn) + if ( p%PlaneSliceSeriesUn(nt) > 0 ) then + call VTK_Series_Close( p%PlaneSliceSeriesUn(nt), errStat, errMsg ) + end if + end do + end if + + ! Feature 3: close any open .vtp.series sidecars + if (allocated(p%TerrainSliceSeriesUn)) then + do nt = 1, size(p%TerrainSliceSeriesUn) + if ( p%TerrainSliceSeriesUn(nt) > 0 ) then + call VTK_Series_Close( p%TerrainSliceSeriesUn(nt), errStat, errMsg ) + end if + end do + end if + ! Destroy InflowWind data select case(p%Mod_AmbWind) case (2) @@ -2004,7 +2980,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], & @@ -2021,7 +2997,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)], & @@ -2038,7 +3014,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)], & @@ -2046,10 +3022,39 @@ subroutine AWAE_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg p%LowRes%dXYZ, m%outVizXZPlane, ErrStat2, ErrMsg2) 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 + + end if + + ! Feature 2: axis-aligned planar sampling (extent-controlled, own sampling rate) + if ( p%NumPlaneSlices > 0 .and. mod(n, p%WrPlaneSkp) == 0 ) then + call EmitPlaneSlices( t, n, p, m, ErrStat2, ErrMsg2 ) + if (Failed()) return + end if + + ! Feature 3: terrain-following point-cloud sampling (own sampling rate) + if ( p%NumTerrainSlices > 0 .and. mod(n, p%WrTerrainSkp) == 0 ) then + call EmitTerrainSlices( t, n, p, m, ErrStat2, ErrMsg2 ) + if (Failed()) return end if contains + logical function Failed() call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) Failed = ErrStat >= AbortErrLev @@ -2581,5 +3586,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..e46bc25a48 100644 --- a/modules/awae/src/AWAE_IO.f90 +++ b/modules/awae/src/AWAE_IO.f90 @@ -36,6 +36,9 @@ MODULE AWAE_IO public :: AWAE_IO_InitGridInfo public :: ReadLowResWindVTK, ReadWindAMReX + public :: WriteVTK_PolyData + public :: WriteVTK_StructuredGrid_2D + public :: VTK_Series_Open, VTK_Series_Append, VTK_Series_Close interface subroutine ReadVTK_inflow_info(FileName, Desc, dims, origin, gridSpacing, vecLabel, values, read_values, err_stat, err_msg) BIND(C,name='ReadVTK_inflow_info') @@ -80,7 +83,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 +95,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 @@ -106,6 +109,233 @@ subroutine WriteDisWindFiles( n, WrDisSkp1, p, y, m, errStat, errMsg ) end subroutine WriteDisWindFiles +!---------------------------------------------------------------------------------------------------------------------------------- +!> Write a single-piece VTK XML `PolyData` (.vtp) file containing an unstructured +!! point cloud with one vertex cell per point plus a 3-component `Vec` field on +!! the points. NaN values in `Vec` are written as the literal text "nan" (VTK +!! ASCII readers accept this and ParaView masks it automatically). +!! +!! Used by `NumTerrainSlices` (Feature 1, terrain-following point cloud). +subroutine WriteVTK_PolyData( FileName, descr, Pts, Vec, vecLabel, ErrStat, ErrMsg ) + character(*), intent(in ) :: FileName !< output filename (usually "...vtp") + character(*), intent(in ) :: descr !< short descriptor written as an XML comment + real(ReKi), intent(in ) :: Pts(:,:) !< (3, N) point coordinates in the FAST.Farm global frame + real(SiKi), intent(in ) :: Vec(:,:) !< (3, N) vector values at each point; NaN is permitted + character(*), intent(in ) :: vecLabel !< name of the vector array in the VTK file + integer(IntKi),intent( out) :: ErrStat !< error status + character(*), intent( out) :: ErrMsg !< error message + + integer(IntKi) :: Un, ErrStat2, i, N + character(ErrMsgLen) :: ErrMsg2 + character(*), parameter :: RoutineName = 'WriteVTK_PolyData' + + ErrStat = ErrID_None + ErrMsg = '' + + N = size(Pts,2) + if (N <= 0) return + if (size(Vec,2) /= N) then + call SetErrStat(ErrID_Fatal, 'Pts and Vec must have matching second dimension', ErrStat, ErrMsg, RoutineName) + return + end if + + !$OMP critical(fileopenNWTCio_critical) + call GetNewUnit( Un, ErrStat2, ErrMsg2 ) + call OpenFOutFile( Un, trim(FileName), ErrStat2, ErrMsg2 ) + !$OMP end critical(fileopenNWTCio_critical) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + + write(Un,'(A)') '' + write(Un,'(A)') '' + write(Un,'(A)') '' + write(Un,'(A)') ' ' + write(Un,'(A,I0,A,I0,A)') ' ' + + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + do i = 1, N + write(Un,'(3(1X,ES14.6))') real(Pts(1,i),SiKi), real(Pts(2,i),SiKi), real(Pts(3,i),SiKi) + end do + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + write(Un,'(*(1X,I0))') (i-1, i=1,N) + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + write(Un,'(*(1X,I0))') (i, i=1,N) + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + do i = 1, N + write(Un,'(3(1X,ES14.6))') Vec(1,i), Vec(2,i), Vec(3,i) + end do + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + write(Un,'(A)') '' + + !$OMP critical(fileopenNWTCio_critical) + close(Un) + !$OMP end critical(fileopenNWTCio_critical) + +end subroutine WriteVTK_PolyData +!---------------------------------------------------------------------------------------------------------------------------------- +!> Write a 2-D `StructuredGrid` (.vts) VTK XML file: an n1 x n2 x 1 mesh with +!! explicit point coordinates plus a 3-component `Vec` field. NaN values are +!! written as the literal text "nan"; ParaView masks these automatically. +!! Node ordering is (i-fastest, then j) matching VTK's Fortran-friendly linear +!! layout. +!! +!! Used by `NumPlaneSlices` (Feature 2, axis-aligned planar sampling) when +!! `VTKFormat=xml`. +subroutine WriteVTK_StructuredGrid_2D( FileName, descr, n1, n2, Pts, Vec, vecLabel, ErrStat, ErrMsg ) + character(*), intent(in ) :: FileName !< output filename + character(*), intent(in ) :: descr !< short descriptor written as an XML comment + integer(IntKi),intent(in ) :: n1, n2 !< dimensions of the 2-D grid + real(ReKi), intent(in ) :: Pts(:,:,:) !< (3, n1, n2) point coordinates in the global frame + real(SiKi), intent(in ) :: Vec(:,:,:) !< (3, n1, n2) vector values at each node; NaN is permitted + character(*), intent(in ) :: vecLabel !< name of the vector array + integer(IntKi),intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + + integer(IntKi) :: Un, ErrStat2, i, j + character(ErrMsgLen) :: ErrMsg2 + character(*), parameter :: RoutineName = 'WriteVTK_StructuredGrid_2D' + character(64) :: extentStr + + ErrStat = ErrID_None + ErrMsg = '' + + if (n1 < 1 .or. n2 < 1) return + + !$OMP critical(fileopenNWTCio_critical) + call GetNewUnit( Un, ErrStat2, ErrMsg2 ) + call OpenFOutFile( Un, trim(FileName), ErrStat2, ErrMsg2 ) + !$OMP end critical(fileopenNWTCio_critical) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + + write(extentStr,'(A,I0,A,I0,A)') '0 ', n1-1, ' 0 ', n2-1, ' 0 0' + + write(Un,'(A)') '' + write(Un,'(A)') '' + write(Un,'(A)') '' + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + do j = 1, n2 + do i = 1, n1 + write(Un,'(3(1X,ES14.6))') real(Pts(1,i,j),SiKi), real(Pts(2,i,j),SiKi), real(Pts(3,i,j),SiKi) + end do + end do + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + do j = 1, n2 + do i = 1, n1 + write(Un,'(3(1X,ES14.6))') Vec(1,i,j), Vec(2,i,j), Vec(3,i,j) + end do + end do + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + write(Un,'(A)') '' + + !$OMP critical(fileopenNWTCio_critical) + close(Un) + !$OMP end critical(fileopenNWTCio_critical) + +end subroutine WriteVTK_StructuredGrid_2D +!---------------------------------------------------------------------------------------------------------------------------------- +!> Open a ParaView `.series` JSON sidecar and write its header. The caller is +!! responsible for appending one entry per output step via `VTK_Series_Append` +!! and closing the file via `VTK_Series_Close`. Used by both `NumPlaneSlices` +!! and `NumTerrainSlices`. +subroutine VTK_Series_Open( FileName, Un, ErrStat, ErrMsg ) + character(*), intent(in ) :: FileName !< output filename (usually ends in ".vts.series" or ".vtp.series") + integer(IntKi),intent( out) :: Un !< Fortran unit number (must be passed unchanged to Append/Close) + integer(IntKi),intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + character(*), parameter :: RoutineName = 'VTK_Series_Open' + + !$OMP critical(fileopenNWTCio_critical) + call GetNewUnit( Un, ErrStat2, ErrMsg2 ) + call OpenFOutFile( Un, trim(FileName), ErrStat2, ErrMsg2 ) + !$OMP end critical(fileopenNWTCio_critical) + ErrStat = ErrID_None; ErrMsg = '' + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + + write(Un,'(A)') '{' + write(Un,'(A)') ' "file-series-version" : "1.0",' + write(Un,'(A)') ' "files" : [' + +end subroutine VTK_Series_Open +!---------------------------------------------------------------------------------------------------------------------------------- +!> Append one file entry to an open `.series` sidecar. `isFirst` must be +!! `.true.` on the very first call and `.false.` thereafter so that the JSON +!! comma placement is valid. +subroutine VTK_Series_Append( Un, fileName, t, isFirst, ErrStat, ErrMsg ) + integer(IntKi),intent(in ) :: Un !< unit returned by VTK_Series_Open + character(*), intent(in ) :: fileName !< name of the VTK file to record (without path) + real(DbKi), intent(in ) :: t !< time in seconds + logical, intent(in ) :: isFirst !< true on the first append, false otherwise + integer(IntKi),intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + + character(:), allocatable :: prefix + + ErrStat = ErrID_None + ErrMsg = '' + + if (isFirst) then + prefix = ' ' + else + prefix = ' ,' + end if + + write(Un,'(A,A,ES14.6,A)') prefix, '{ "name" : "'//trim(fileName)//'", "time" : ', t, ' }' + +end subroutine VTK_Series_Append +!---------------------------------------------------------------------------------------------------------------------------------- +!> Close an open `.series` sidecar. Writes the terminating `]}` and closes +!! the Fortran unit. +subroutine VTK_Series_Close( Un, ErrStat, ErrMsg ) + integer(IntKi),intent(inout) :: Un !< unit returned by VTK_Series_Open; set to -1 on exit + integer(IntKi),intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + + ErrStat = ErrID_None + ErrMsg = '' + + if (Un <= 0) return + + write(Un,'(A)') ' ]' + write(Un,'(A)') '}' + + !$OMP critical(fileopenNWTCio_critical) + close(Un) + !$OMP end critical(fileopenNWTCio_critical) + Un = -1 + +end subroutine VTK_Series_Close + !---------------------------------------------------------------------------------------------------------------------------------- !> This subroutine read the low res wind file (VTK) at a given time step `n` subroutine ReadLowResWindVTK(n, p, Vamb_Low, errStat, errMsg) @@ -328,7 +558,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 diff --git a/modules/awae/src/AWAE_Registry.txt b/modules/awae/src/AWAE_Registry.txt index 7a52d7e940..425d020a47 100644 --- a/modules/awae/src/AWAE_Registry.txt +++ b/modules/awae/src/AWAE_Registry.txt @@ -40,6 +40,38 @@ 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 + +# ..... Axis-aligned planar sampling (Feature 2, Snippet 6 grammar) ....................................................... +# Feature 2 uses flat arrays (rather than a derived slice type) so that the +# generated pack/unpack routines stay small enough to compile under 8 GB of +# container RAM. Populated by the FAST.Farm input reader from the +# '--- AXIS-ALIGNED PLANE SLICES (extent-controlled) ---' block, then consumed +# by AWAE_Init to compute derived indices and by AWAE_CalcOutput to emit +# per-step VTK output. See modules/awae/src/AWAE.f90 for the semantics of the +# derived NormalAxis / N1 / N2 / bounds fields. +typedef ^ AWAE_InputFileType DbKi WrPlaneDT - - - "Feature 2 sampling period (s); DEFAULT falls back to DT_Low" s +typedef ^ AWAE_InputFileType IntKi NumPlaneSlices - - - "Number of axis-aligned extent-controlled planar slices (Feature 2)" - +typedef ^ AWAE_InputFileType CHARACTER(64) PlaneSliceName {:} - - "user-provided slice name (used in output filenames)" - +typedef ^ AWAE_InputFileType ReKi PlaneSliceOrigin {3}{:} - - "plane corner (x0,y0,z0) in the global farm frame; column k is slice k" m +typedef ^ AWAE_InputFileType ReKi PlaneSliceNormal {3}{:} - - "plane normal derived from Plane keyword (XY→001, YZ→100, XZ→010); column k is slice k" - +typedef ^ AWAE_InputFileType ReKi PlaneSliceExtent1 {:} - - "in-plane extent along the first non-normal axis, per slice" m +typedef ^ AWAE_InputFileType ReKi PlaneSliceExtent2 {:} - - "in-plane extent along the second non-normal axis, per slice" m + +# ..... Terrain-following point-cloud sampling (Feature 3) ....................................................... +# Feature 3 stores the raw config (name, source flag, filename, offsets, normal) in the input file type. +# Offsets are ragged (each slice can have a different number of offsets) but stored via a flat/packed +# layout: OffsetsPacked contains all offsets concatenated, OffsetIdx(k) marks the start of slice k in +# the packed array with the convention that slice k occupies indices [OffsetIdx(k)+1 .. OffsetIdx(k+1)]. +# See modules/awae/src/AWAE.f90 for the parser and the point-cloud/STL readers. +typedef ^ AWAE_InputFileType DbKi WrTerrainDT - - - "Feature 3 sampling period (s); DEFAULT falls back to DT_Low" s +typedef ^ AWAE_InputFileType IntKi NumTerrainSlices - - - "Number of terrain-following point-cloud slices (Feature 3)" - +typedef ^ AWAE_InputFileType CHARACTER(64) TerrainSliceName {:} - - "user-provided slice name (used in output filenames)" - +typedef ^ AWAE_InputFileType IntKi TerrainSliceSourceType {:} - - "1 = STL, 2 = plain-text point cloud" - +typedef ^ AWAE_InputFileType CHARACTER(1024) TerrainSliceFileName {:} - - "STL or point-cloud file name; resolved relative to the primary input file" - +typedef ^ AWAE_InputFileType ReKi TerrainSliceOffsetNormal {3}{:} - - "unit vector along which each slice's Offsets are applied; if all zero, per-facet normal from STL is used" - +typedef ^ AWAE_InputFileType ReKi TerrainSliceOffsetsPacked {:} - - "packed per-slice offsets (m); slice k occupies TerrainSliceOffsetIdx(k)+1 .. TerrainSliceOffsetIdx(k+1)" m +typedef ^ AWAE_InputFileType IntKi TerrainSliceOffsetIdx {:} - - "prefix-sum indices into OffsetsPacked; size = NumTerrainSlices+1" - + 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,15 +159,19 @@ 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 IntKi PlaneSliceStepCount {:} - - "Feature 2: number of steps written per slice (drives isFirst flag in the .vts.series writer)" - +typedef ^ MiscVarType IntKi TerrainSliceStepCount {:} - - "Feature 3: number of steps written per slice (drives isFirst flag in the .vtp.series writer)" - +typedef ^ MiscVarType IntKi LastPlaneSliceN - -1 - "Feature 2: last low-res step index at which planar slices were emitted; -1 initially. Guards the duplicate-emit at t=0 triggered by FARM_InitialCO calling AWAE_CalcOutput twice." - +typedef ^ MiscVarType IntKi LastTerrainSliceN - -1 - "Feature 3: last low-res step index at which terrain slices were emitted; -1 initially. Same guard purpose as LastPlaneSliceN." - typedef ^ MiscVarType SiKi outVizYZPlane {:}{:}{:}{:} - - "An array holding the output data for a 2D visualization slice" - typedef ^ MiscVarType SiKi outVizXZPlane {:}{:}{:}{:} - - "An array holding the output data for a 2D visualization slice" - typedef ^ MiscVarType InflowWind_InputType u_IfW_Low - - - "InflowWind module inputs for the low-resolution grid" - @@ -145,7 +181,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 {:}{:} - - "Time step when wake plane starts - counted by N_dtLow. Indices [wakenum,turbnum]" - # Low-resolution grid chunk data typedef ^ LRGChunkType IntKi iChunk {3} - - "XYZ index of chunk" - @@ -166,7 +205,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" - @@ -224,9 +262,38 @@ typedef ^ ParameterType LOGICAL OutDisWindXvalid {:} - - "Valid Y typedef ^ ParameterType 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 ^ 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]" - +# Feature 2 runtime state (flat arrays, indexed by slice number 1..NumPlaneSlices). +typedef ^ ParameterType IntKi WrPlaneSkp - - - "Feature 2: number of low-res time steps between planar sampling emits (WrPlaneDT/dt_low)" - +typedef ^ ParameterType IntKi NumPlaneSlices - - - "Feature 2: number of axis-aligned extent-controlled slices" - +typedef ^ ParameterType CHARACTER(64) PlaneSliceName {:} - - "Feature 2: per-slice user name" - +typedef ^ ParameterType ReKi PlaneSliceOrigin {3}{:} - - "Feature 2: (3, NumPlaneSlices) plane corner in the global farm frame" m +typedef ^ ParameterType IntKi PlaneSliceNormalAxis {:} - - "Feature 2: 1=X, 2=Y, 3=Z; derived from Normal at init" - +typedef ^ ParameterType IntKi PlaneSliceAxis1 {:} - - "Feature 2: index (1..3) of first in-plane global axis" - +typedef ^ ParameterType IntKi PlaneSliceAxis2 {:} - - "Feature 2: index (1..3) of second in-plane global axis" - +typedef ^ ParameterType IntKi PlaneSliceN1 {:} - - "Feature 2: number of in-plane nodes along Axis1" - +typedef ^ ParameterType IntKi PlaneSliceN2 {:} - - "Feature 2: number of in-plane nodes along Axis2" - +typedef ^ ParameterType IntKi PlaneSliceILoReq {:} - - "Feature 2: 0-based lo-index along Axis1 in parent low-res grid (may be negative)" - +typedef ^ ParameterType IntKi PlaneSliceIHiReq {:} - - "Feature 2: 0-based hi-index along Axis1 (may exceed nAxis1-1)" - +typedef ^ ParameterType IntKi PlaneSliceJLoReq {:} - - "Feature 2: 0-based lo-index along Axis2" - +typedef ^ ParameterType IntKi PlaneSliceJHiReq {:} - - "Feature 2: 0-based hi-index along Axis2" - +typedef ^ ParameterType LOGICAL PlaneSliceValid {:} - - "Feature 2: true if the slice intersects the low-res domain" - +typedef ^ ParameterType IntKi PlaneSliceSeriesUn {:} - - "Feature 2: open .vts.series unit per slice; -1 = closed" - +# Feature 3 runtime state (flat arrays, indexed by slice number 1..NumTerrainSlices). +typedef ^ ParameterType IntKi WrTerrainSkp - - - "Feature 3: number of low-res time steps between terrain sampling emits (WrTerrainDT/dt_low)" - +typedef ^ ParameterType IntKi NumTerrainSlices - - - "Feature 3: number of terrain-following point-cloud slices" - +typedef ^ ParameterType CHARACTER(64) TerrainSliceName {:} - - "Feature 3: per-slice user name" - +typedef ^ ParameterType LOGICAL TerrainSliceValid {:} - - "Feature 3: true if the slice was read successfully and has at least one in-domain point" - +typedef ^ ParameterType IntKi TerrainSliceSeriesUn {:} - - "Feature 3: open .vtp.series unit per slice; -1 = closed" - +typedef ^ ParameterType IntKi TerrainSliceNPtsTotal {:} - - "Feature 3: total number of output points per slice (base NPts x N_offsets)" - +typedef ^ ParameterType IntKi TerrainSlicePtsIdx {:} - - "Feature 3: prefix-sum indices into TerrainSlicePtsFlat; size = NumTerrainSlices+1" - +typedef ^ ParameterType ReKi TerrainSlicePtsFlat {3}{:} - - "Feature 3: concatenated per-slice point coordinates in the farm-global frame; slice k occupies TerrainSlicePtsIdx(k)+1 .. TerrainSlicePtsIdx(k+1)" m 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..eb2c11d4be 100644 --- a/modules/awae/src/AWAE_Types.f90 +++ b/modules/awae/src/AWAE_Types.f90 @@ -65,6 +65,21 @@ 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] + REAL(DbKi) :: WrPlaneDT = 0.0_R8Ki !< Feature 2 sampling period (s); DEFAULT falls back to DT_Low [s] + INTEGER(IntKi) :: NumPlaneSlices = 0_IntKi !< Number of axis-aligned extent-controlled planar slices (Feature 2) [-] + CHARACTER(64) , DIMENSION(:), ALLOCATABLE :: PlaneSliceName !< user-provided slice name (used in output filenames) [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: PlaneSliceOrigin !< plane corner (x0,y0,z0) in the global farm frame; column k is slice k [m] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: PlaneSliceNormal !< plane normal derived from Plane keyword (XY→001, YZ→100, XZ→010); column k is slice k [-] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceExtent1 !< in-plane extent along the first non-normal axis, per slice [m] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceExtent2 !< in-plane extent along the second non-normal axis, per slice [m] + REAL(DbKi) :: WrTerrainDT = 0.0_R8Ki !< Feature 3 sampling period (s); DEFAULT falls back to DT_Low [s] + INTEGER(IntKi) :: NumTerrainSlices = 0_IntKi !< Number of terrain-following point-cloud slices (Feature 3) [-] + CHARACTER(64) , DIMENSION(:), ALLOCATABLE :: TerrainSliceName !< user-provided slice name (used in output filenames) [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: TerrainSliceSourceType !< 1 = STL, 2 = plain-text point cloud [-] + CHARACTER(1024) , DIMENSION(:), ALLOCATABLE :: TerrainSliceFileName !< STL or point-cloud file name; resolved relative to the primary input file [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: TerrainSliceOffsetNormal !< unit vector along which each slice's Offsets are applied; if all zero, per-facet normal from STL is used [-] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: TerrainSliceOffsetsPacked !< packed per-slice offsets (m); slice k occupies TerrainSliceOffsetIdx(k)+1 .. TerrainSliceOffsetIdx(k+1) [m] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: TerrainSliceOffsetIdx !< prefix-sum indices into OffsetsPacked; size = NumTerrainSlices+1 [-] 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,14 +169,18 @@ 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 [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceStepCount !< Feature 2: number of steps written per slice (drives isFirst flag in the .vts.series writer) [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: TerrainSliceStepCount !< Feature 3: number of steps written per slice (drives isFirst flag in the .vtp.series writer) [-] + INTEGER(IntKi) :: LastPlaneSliceN = -1 !< Feature 2: last low-res step index at which planar slices were emitted; -1 initially. Guards the duplicate-emit at t=0 triggered by FARM_InitialCO calling AWAE_CalcOutput twice. [-] + INTEGER(IntKi) :: LastTerrainSliceN = -1 !< Feature 3: last low-res step index at which terrain slices were emitted; -1 initially. Same guard purpose as LastPlaneSliceN. [-] 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 [-] TYPE(InflowWind_InputType) :: u_IfW_Low !< InflowWind module inputs for the low-resolution grid [-] @@ -169,7 +188,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 !< Time step when wake plane starts - counted by N_dtLow. Indices [wakenum,turbnum] [-] END TYPE AWAE_MiscVarType ! ======================= ! ========= LRGChunkType ======= @@ -194,7 +214,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 ! ======================= @@ -250,9 +269,36 @@ 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] LOGICAL , DIMENSION(:), ALLOCATABLE :: OutDisWindYvalid !< Valid XZ planes for output of disturbed wind data across the low-resolution domain [1 to NOutDisWindXZ] [-] + INTEGER(IntKi) :: WrPlaneSkp = 0_IntKi !< Feature 2: number of low-res time steps between planar sampling emits (WrPlaneDT/dt_low) [-] + INTEGER(IntKi) :: NumPlaneSlices = 0_IntKi !< Feature 2: number of axis-aligned extent-controlled slices [-] + CHARACTER(64) , DIMENSION(:), ALLOCATABLE :: PlaneSliceName !< Feature 2: per-slice user name [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: PlaneSliceOrigin !< Feature 2: (3, NumPlaneSlices) plane corner in the global farm frame [m] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceNormalAxis !< Feature 2: 1=X, 2=Y, 3=Z; derived from Normal at init [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceAxis1 !< Feature 2: index (1..3) of first in-plane global axis [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceAxis2 !< Feature 2: index (1..3) of second in-plane global axis [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceN1 !< Feature 2: number of in-plane nodes along Axis1 [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceN2 !< Feature 2: number of in-plane nodes along Axis2 [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceILoReq !< Feature 2: 0-based lo-index along Axis1 in parent low-res grid (may be negative) [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceIHiReq !< Feature 2: 0-based hi-index along Axis1 (may exceed nAxis1-1) [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceJLoReq !< Feature 2: 0-based lo-index along Axis2 [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceJHiReq !< Feature 2: 0-based hi-index along Axis2 [-] + LOGICAL , DIMENSION(:), ALLOCATABLE :: PlaneSliceValid !< Feature 2: true if the slice intersects the low-res domain [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: PlaneSliceSeriesUn !< Feature 2: open .vts.series unit per slice; -1 = closed [-] + INTEGER(IntKi) :: WrTerrainSkp = 0_IntKi !< Feature 3: number of low-res time steps between terrain sampling emits (WrTerrainDT/dt_low) [-] + INTEGER(IntKi) :: NumTerrainSlices = 0_IntKi !< Feature 3: number of terrain-following point-cloud slices [-] + CHARACTER(64) , DIMENSION(:), ALLOCATABLE :: TerrainSliceName !< Feature 3: per-slice user name [-] + LOGICAL , DIMENSION(:), ALLOCATABLE :: TerrainSliceValid !< Feature 3: true if the slice was read successfully and has at least one in-domain point [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: TerrainSliceSeriesUn !< Feature 3: open .vtp.series unit per slice; -1 = closed [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: TerrainSliceNPtsTotal !< Feature 3: total number of output points per slice (base NPts x N_offsets) [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: TerrainSlicePtsIdx !< Feature 3: prefix-sum indices into TerrainSlicePtsFlat; size = NumTerrainSlices+1 [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: TerrainSlicePtsFlat !< Feature 3: concatenated per-slice point coordinates in the farm-global frame; slice k occupies TerrainSlicePtsIdx(k)+1 .. TerrainSlicePtsIdx(k+1) [m] 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 +504,142 @@ subroutine AWAE_CopyInputFileType(SrcInputFileTypeData, DstInputFileTypeData, Ct DstInputFileTypeData%OutDisWindY = SrcInputFileTypeData%OutDisWindY end if DstInputFileTypeData%WrDisDT = SrcInputFileTypeData%WrDisDT + DstInputFileTypeData%WrPlaneDT = SrcInputFileTypeData%WrPlaneDT + DstInputFileTypeData%NumPlaneSlices = SrcInputFileTypeData%NumPlaneSlices + if (allocated(SrcInputFileTypeData%PlaneSliceName)) then + LB(1:1) = lbound(SrcInputFileTypeData%PlaneSliceName) + UB(1:1) = ubound(SrcInputFileTypeData%PlaneSliceName) + if (.not. allocated(DstInputFileTypeData%PlaneSliceName)) then + allocate(DstInputFileTypeData%PlaneSliceName(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstInputFileTypeData%PlaneSliceName.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstInputFileTypeData%PlaneSliceName = SrcInputFileTypeData%PlaneSliceName + end if + if (allocated(SrcInputFileTypeData%PlaneSliceOrigin)) then + LB(1:2) = lbound(SrcInputFileTypeData%PlaneSliceOrigin) + UB(1:2) = ubound(SrcInputFileTypeData%PlaneSliceOrigin) + if (.not. allocated(DstInputFileTypeData%PlaneSliceOrigin)) then + allocate(DstInputFileTypeData%PlaneSliceOrigin(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstInputFileTypeData%PlaneSliceOrigin.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstInputFileTypeData%PlaneSliceOrigin = SrcInputFileTypeData%PlaneSliceOrigin + end if + if (allocated(SrcInputFileTypeData%PlaneSliceNormal)) then + LB(1:2) = lbound(SrcInputFileTypeData%PlaneSliceNormal) + UB(1:2) = ubound(SrcInputFileTypeData%PlaneSliceNormal) + if (.not. allocated(DstInputFileTypeData%PlaneSliceNormal)) then + allocate(DstInputFileTypeData%PlaneSliceNormal(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstInputFileTypeData%PlaneSliceNormal.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstInputFileTypeData%PlaneSliceNormal = SrcInputFileTypeData%PlaneSliceNormal + end if + if (allocated(SrcInputFileTypeData%PlaneSliceExtent1)) then + LB(1:1) = lbound(SrcInputFileTypeData%PlaneSliceExtent1) + UB(1:1) = ubound(SrcInputFileTypeData%PlaneSliceExtent1) + if (.not. allocated(DstInputFileTypeData%PlaneSliceExtent1)) then + allocate(DstInputFileTypeData%PlaneSliceExtent1(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstInputFileTypeData%PlaneSliceExtent1.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstInputFileTypeData%PlaneSliceExtent1 = SrcInputFileTypeData%PlaneSliceExtent1 + end if + if (allocated(SrcInputFileTypeData%PlaneSliceExtent2)) then + LB(1:1) = lbound(SrcInputFileTypeData%PlaneSliceExtent2) + UB(1:1) = ubound(SrcInputFileTypeData%PlaneSliceExtent2) + if (.not. allocated(DstInputFileTypeData%PlaneSliceExtent2)) then + allocate(DstInputFileTypeData%PlaneSliceExtent2(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstInputFileTypeData%PlaneSliceExtent2.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstInputFileTypeData%PlaneSliceExtent2 = SrcInputFileTypeData%PlaneSliceExtent2 + end if + DstInputFileTypeData%WrTerrainDT = SrcInputFileTypeData%WrTerrainDT + DstInputFileTypeData%NumTerrainSlices = SrcInputFileTypeData%NumTerrainSlices + if (allocated(SrcInputFileTypeData%TerrainSliceName)) then + LB(1:1) = lbound(SrcInputFileTypeData%TerrainSliceName) + UB(1:1) = ubound(SrcInputFileTypeData%TerrainSliceName) + if (.not. allocated(DstInputFileTypeData%TerrainSliceName)) then + allocate(DstInputFileTypeData%TerrainSliceName(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstInputFileTypeData%TerrainSliceName.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstInputFileTypeData%TerrainSliceName = SrcInputFileTypeData%TerrainSliceName + end if + if (allocated(SrcInputFileTypeData%TerrainSliceSourceType)) then + LB(1:1) = lbound(SrcInputFileTypeData%TerrainSliceSourceType) + UB(1:1) = ubound(SrcInputFileTypeData%TerrainSliceSourceType) + if (.not. allocated(DstInputFileTypeData%TerrainSliceSourceType)) then + allocate(DstInputFileTypeData%TerrainSliceSourceType(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstInputFileTypeData%TerrainSliceSourceType.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstInputFileTypeData%TerrainSliceSourceType = SrcInputFileTypeData%TerrainSliceSourceType + end if + if (allocated(SrcInputFileTypeData%TerrainSliceFileName)) then + LB(1:1) = lbound(SrcInputFileTypeData%TerrainSliceFileName) + UB(1:1) = ubound(SrcInputFileTypeData%TerrainSliceFileName) + if (.not. allocated(DstInputFileTypeData%TerrainSliceFileName)) then + allocate(DstInputFileTypeData%TerrainSliceFileName(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstInputFileTypeData%TerrainSliceFileName.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstInputFileTypeData%TerrainSliceFileName = SrcInputFileTypeData%TerrainSliceFileName + end if + if (allocated(SrcInputFileTypeData%TerrainSliceOffsetNormal)) then + LB(1:2) = lbound(SrcInputFileTypeData%TerrainSliceOffsetNormal) + UB(1:2) = ubound(SrcInputFileTypeData%TerrainSliceOffsetNormal) + if (.not. allocated(DstInputFileTypeData%TerrainSliceOffsetNormal)) then + allocate(DstInputFileTypeData%TerrainSliceOffsetNormal(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstInputFileTypeData%TerrainSliceOffsetNormal.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstInputFileTypeData%TerrainSliceOffsetNormal = SrcInputFileTypeData%TerrainSliceOffsetNormal + end if + if (allocated(SrcInputFileTypeData%TerrainSliceOffsetsPacked)) then + LB(1:1) = lbound(SrcInputFileTypeData%TerrainSliceOffsetsPacked) + UB(1:1) = ubound(SrcInputFileTypeData%TerrainSliceOffsetsPacked) + if (.not. allocated(DstInputFileTypeData%TerrainSliceOffsetsPacked)) then + allocate(DstInputFileTypeData%TerrainSliceOffsetsPacked(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstInputFileTypeData%TerrainSliceOffsetsPacked.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstInputFileTypeData%TerrainSliceOffsetsPacked = SrcInputFileTypeData%TerrainSliceOffsetsPacked + end if + if (allocated(SrcInputFileTypeData%TerrainSliceOffsetIdx)) then + LB(1:1) = lbound(SrcInputFileTypeData%TerrainSliceOffsetIdx) + UB(1:1) = ubound(SrcInputFileTypeData%TerrainSliceOffsetIdx) + if (.not. allocated(DstInputFileTypeData%TerrainSliceOffsetIdx)) then + allocate(DstInputFileTypeData%TerrainSliceOffsetIdx(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstInputFileTypeData%TerrainSliceOffsetIdx.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstInputFileTypeData%TerrainSliceOffsetIdx = SrcInputFileTypeData%TerrainSliceOffsetIdx + end if DstInputFileTypeData%ChkWndFiles = SrcInputFileTypeData%ChkWndFiles DstInputFileTypeData%Mod_Meander = SrcInputFileTypeData%Mod_Meander DstInputFileTypeData%C_Meander = SrcInputFileTypeData%C_Meander @@ -580,6 +762,39 @@ subroutine AWAE_DestroyInputFileType(InputFileTypeData, ErrStat, ErrMsg) if (allocated(InputFileTypeData%OutDisWindY)) then deallocate(InputFileTypeData%OutDisWindY) end if + if (allocated(InputFileTypeData%PlaneSliceName)) then + deallocate(InputFileTypeData%PlaneSliceName) + end if + if (allocated(InputFileTypeData%PlaneSliceOrigin)) then + deallocate(InputFileTypeData%PlaneSliceOrigin) + end if + if (allocated(InputFileTypeData%PlaneSliceNormal)) then + deallocate(InputFileTypeData%PlaneSliceNormal) + end if + if (allocated(InputFileTypeData%PlaneSliceExtent1)) then + deallocate(InputFileTypeData%PlaneSliceExtent1) + end if + if (allocated(InputFileTypeData%PlaneSliceExtent2)) then + deallocate(InputFileTypeData%PlaneSliceExtent2) + end if + if (allocated(InputFileTypeData%TerrainSliceName)) then + deallocate(InputFileTypeData%TerrainSliceName) + end if + if (allocated(InputFileTypeData%TerrainSliceSourceType)) then + deallocate(InputFileTypeData%TerrainSliceSourceType) + end if + if (allocated(InputFileTypeData%TerrainSliceFileName)) then + deallocate(InputFileTypeData%TerrainSliceFileName) + end if + if (allocated(InputFileTypeData%TerrainSliceOffsetNormal)) then + deallocate(InputFileTypeData%TerrainSliceOffsetNormal) + end if + if (allocated(InputFileTypeData%TerrainSliceOffsetsPacked)) then + deallocate(InputFileTypeData%TerrainSliceOffsetsPacked) + end if + if (allocated(InputFileTypeData%TerrainSliceOffsetIdx)) then + deallocate(InputFileTypeData%TerrainSliceOffsetIdx) + end if if (allocated(InputFileTypeData%X0_high)) then deallocate(InputFileTypeData%X0_high) end if @@ -621,6 +836,21 @@ subroutine AWAE_PackInputFileType(RF, Indata) call RegPack(RF, InData%NOutDisWindXZ) call RegPackAlloc(RF, InData%OutDisWindY) call RegPack(RF, InData%WrDisDT) + call RegPack(RF, InData%WrPlaneDT) + call RegPack(RF, InData%NumPlaneSlices) + call RegPackAlloc(RF, InData%PlaneSliceName) + call RegPackAlloc(RF, InData%PlaneSliceOrigin) + call RegPackAlloc(RF, InData%PlaneSliceNormal) + call RegPackAlloc(RF, InData%PlaneSliceExtent1) + call RegPackAlloc(RF, InData%PlaneSliceExtent2) + call RegPack(RF, InData%WrTerrainDT) + call RegPack(RF, InData%NumTerrainSlices) + call RegPackAlloc(RF, InData%TerrainSliceName) + call RegPackAlloc(RF, InData%TerrainSliceSourceType) + call RegPackAlloc(RF, InData%TerrainSliceFileName) + call RegPackAlloc(RF, InData%TerrainSliceOffsetNormal) + call RegPackAlloc(RF, InData%TerrainSliceOffsetsPacked) + call RegPackAlloc(RF, InData%TerrainSliceOffsetIdx) call RegPack(RF, InData%ChkWndFiles) call RegPack(RF, InData%Mod_Meander) call RegPack(RF, InData%C_Meander) @@ -672,6 +902,21 @@ 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%WrPlaneDT); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%NumPlaneSlices); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceName); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceOrigin); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceNormal); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceExtent1); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceExtent2); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%WrTerrainDT); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%NumTerrainSlices); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSliceName); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSliceSourceType); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSliceFileName); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSliceOffsetNormal); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSliceOffsetsPacked); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSliceOffsetIdx); 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 @@ -1366,6 +1611,32 @@ subroutine AWAE_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) end if DstMiscData%outVizXYPlane = SrcMiscData%outVizXYPlane end if + if (allocated(SrcMiscData%PlaneSliceStepCount)) then + LB(1:1) = lbound(SrcMiscData%PlaneSliceStepCount) + UB(1:1) = ubound(SrcMiscData%PlaneSliceStepCount) + if (.not. allocated(DstMiscData%PlaneSliceStepCount)) then + allocate(DstMiscData%PlaneSliceStepCount(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%PlaneSliceStepCount.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstMiscData%PlaneSliceStepCount = SrcMiscData%PlaneSliceStepCount + end if + if (allocated(SrcMiscData%TerrainSliceStepCount)) then + LB(1:1) = lbound(SrcMiscData%TerrainSliceStepCount) + UB(1:1) = ubound(SrcMiscData%TerrainSliceStepCount) + if (.not. allocated(DstMiscData%TerrainSliceStepCount)) then + allocate(DstMiscData%TerrainSliceStepCount(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%TerrainSliceStepCount.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstMiscData%TerrainSliceStepCount = SrcMiscData%TerrainSliceStepCount + end if + DstMiscData%LastPlaneSliceN = SrcMiscData%LastPlaneSliceN + DstMiscData%LastTerrainSliceN = SrcMiscData%LastTerrainSliceN if (allocated(SrcMiscData%outVizYZPlane)) then LB(1:4) = lbound(SrcMiscData%outVizYZPlane) UB(1:4) = ubound(SrcMiscData%outVizYZPlane) @@ -1441,10 +1712,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 +1723,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) @@ -1530,6 +1813,12 @@ subroutine AWAE_DestroyMisc(MiscData, ErrStat, ErrMsg) if (allocated(MiscData%outVizXYPlane)) then deallocate(MiscData%outVizXYPlane) end if + if (allocated(MiscData%PlaneSliceStepCount)) then + deallocate(MiscData%PlaneSliceStepCount) + end if + if (allocated(MiscData%TerrainSliceStepCount)) then + deallocate(MiscData%TerrainSliceStepCount) + end if if (allocated(MiscData%outVizYZPlane)) then deallocate(MiscData%outVizYZPlane) end if @@ -1564,6 +1853,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) @@ -1602,6 +1894,10 @@ subroutine AWAE_PackMisc(RF, Indata) call RegPackAlloc(RF, InData%pvec_cs) call RegPackAlloc(RF, InData%pvec_ce) call RegPackAlloc(RF, InData%outVizXYPlane) + call RegPackAlloc(RF, InData%PlaneSliceStepCount) + call RegPackAlloc(RF, InData%TerrainSliceStepCount) + call RegPack(RF, InData%LastPlaneSliceN) + call RegPack(RF, InData%LastTerrainSliceN) call RegPackAlloc(RF, InData%outVizYZPlane) call RegPackAlloc(RF, InData%outVizXZPlane) call InflowWind_PackInput(RF, InData%u_IfW_Low) @@ -1626,6 +1922,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 @@ -1671,6 +1968,10 @@ subroutine AWAE_UnPackMisc(RF, OutData) call RegUnpackAlloc(RF, OutData%pvec_cs); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%pvec_ce); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%outVizXYPlane); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceStepCount); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSliceStepCount); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%LastPlaneSliceN); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%LastTerrainSliceN); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%outVizYZPlane); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%outVizXZPlane); if (RegCheckErr(RF, RoutineName)) return call InflowWind_UnpackInput(RF, OutData%u_IfW_Low) ! u_IfW_Low @@ -1703,6 +2004,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 +2121,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 +2177,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 +2204,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 @@ -2007,8 +2306,8 @@ subroutine AWAE_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) integer(IntKi), intent(in ) :: CtrlCode integer(IntKi), intent( out) :: ErrStat character(*), intent( out) :: ErrMsg - integer(B4Ki) :: i1 - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: i1, i2 + integer(B4Ki) :: LB(2), UB(2) integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 character(*), parameter :: RoutineName = 'AWAE_CopyParam' @@ -2171,9 +2470,245 @@ subroutine AWAE_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) end if DstParamData%OutDisWindYvalid = SrcParamData%OutDisWindYvalid end if + DstParamData%WrPlaneSkp = SrcParamData%WrPlaneSkp + DstParamData%NumPlaneSlices = SrcParamData%NumPlaneSlices + if (allocated(SrcParamData%PlaneSliceName)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceName) + UB(1:1) = ubound(SrcParamData%PlaneSliceName) + if (.not. allocated(DstParamData%PlaneSliceName)) then + allocate(DstParamData%PlaneSliceName(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceName.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceName = SrcParamData%PlaneSliceName + end if + if (allocated(SrcParamData%PlaneSliceOrigin)) then + LB(1:2) = lbound(SrcParamData%PlaneSliceOrigin) + UB(1:2) = ubound(SrcParamData%PlaneSliceOrigin) + if (.not. allocated(DstParamData%PlaneSliceOrigin)) then + allocate(DstParamData%PlaneSliceOrigin(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceOrigin.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceOrigin = SrcParamData%PlaneSliceOrigin + end if + if (allocated(SrcParamData%PlaneSliceNormalAxis)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceNormalAxis) + UB(1:1) = ubound(SrcParamData%PlaneSliceNormalAxis) + if (.not. allocated(DstParamData%PlaneSliceNormalAxis)) then + allocate(DstParamData%PlaneSliceNormalAxis(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceNormalAxis.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceNormalAxis = SrcParamData%PlaneSliceNormalAxis + end if + if (allocated(SrcParamData%PlaneSliceAxis1)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceAxis1) + UB(1:1) = ubound(SrcParamData%PlaneSliceAxis1) + if (.not. allocated(DstParamData%PlaneSliceAxis1)) then + allocate(DstParamData%PlaneSliceAxis1(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceAxis1.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceAxis1 = SrcParamData%PlaneSliceAxis1 + end if + if (allocated(SrcParamData%PlaneSliceAxis2)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceAxis2) + UB(1:1) = ubound(SrcParamData%PlaneSliceAxis2) + if (.not. allocated(DstParamData%PlaneSliceAxis2)) then + allocate(DstParamData%PlaneSliceAxis2(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceAxis2.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceAxis2 = SrcParamData%PlaneSliceAxis2 + end if + if (allocated(SrcParamData%PlaneSliceN1)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceN1) + UB(1:1) = ubound(SrcParamData%PlaneSliceN1) + if (.not. allocated(DstParamData%PlaneSliceN1)) then + allocate(DstParamData%PlaneSliceN1(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceN1.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceN1 = SrcParamData%PlaneSliceN1 + end if + if (allocated(SrcParamData%PlaneSliceN2)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceN2) + UB(1:1) = ubound(SrcParamData%PlaneSliceN2) + if (.not. allocated(DstParamData%PlaneSliceN2)) then + allocate(DstParamData%PlaneSliceN2(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceN2.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceN2 = SrcParamData%PlaneSliceN2 + end if + if (allocated(SrcParamData%PlaneSliceILoReq)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceILoReq) + UB(1:1) = ubound(SrcParamData%PlaneSliceILoReq) + if (.not. allocated(DstParamData%PlaneSliceILoReq)) then + allocate(DstParamData%PlaneSliceILoReq(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceILoReq.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceILoReq = SrcParamData%PlaneSliceILoReq + end if + if (allocated(SrcParamData%PlaneSliceIHiReq)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceIHiReq) + UB(1:1) = ubound(SrcParamData%PlaneSliceIHiReq) + if (.not. allocated(DstParamData%PlaneSliceIHiReq)) then + allocate(DstParamData%PlaneSliceIHiReq(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceIHiReq.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceIHiReq = SrcParamData%PlaneSliceIHiReq + end if + if (allocated(SrcParamData%PlaneSliceJLoReq)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceJLoReq) + UB(1:1) = ubound(SrcParamData%PlaneSliceJLoReq) + if (.not. allocated(DstParamData%PlaneSliceJLoReq)) then + allocate(DstParamData%PlaneSliceJLoReq(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceJLoReq.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceJLoReq = SrcParamData%PlaneSliceJLoReq + end if + if (allocated(SrcParamData%PlaneSliceJHiReq)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceJHiReq) + UB(1:1) = ubound(SrcParamData%PlaneSliceJHiReq) + if (.not. allocated(DstParamData%PlaneSliceJHiReq)) then + allocate(DstParamData%PlaneSliceJHiReq(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceJHiReq.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceJHiReq = SrcParamData%PlaneSliceJHiReq + end if + if (allocated(SrcParamData%PlaneSliceValid)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceValid) + UB(1:1) = ubound(SrcParamData%PlaneSliceValid) + if (.not. allocated(DstParamData%PlaneSliceValid)) then + allocate(DstParamData%PlaneSliceValid(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceValid.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceValid = SrcParamData%PlaneSliceValid + end if + if (allocated(SrcParamData%PlaneSliceSeriesUn)) then + LB(1:1) = lbound(SrcParamData%PlaneSliceSeriesUn) + UB(1:1) = ubound(SrcParamData%PlaneSliceSeriesUn) + if (.not. allocated(DstParamData%PlaneSliceSeriesUn)) then + allocate(DstParamData%PlaneSliceSeriesUn(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%PlaneSliceSeriesUn.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%PlaneSliceSeriesUn = SrcParamData%PlaneSliceSeriesUn + end if + DstParamData%WrTerrainSkp = SrcParamData%WrTerrainSkp + DstParamData%NumTerrainSlices = SrcParamData%NumTerrainSlices + if (allocated(SrcParamData%TerrainSliceName)) then + LB(1:1) = lbound(SrcParamData%TerrainSliceName) + UB(1:1) = ubound(SrcParamData%TerrainSliceName) + if (.not. allocated(DstParamData%TerrainSliceName)) then + allocate(DstParamData%TerrainSliceName(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%TerrainSliceName.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%TerrainSliceName = SrcParamData%TerrainSliceName + end if + if (allocated(SrcParamData%TerrainSliceValid)) then + LB(1:1) = lbound(SrcParamData%TerrainSliceValid) + UB(1:1) = ubound(SrcParamData%TerrainSliceValid) + if (.not. allocated(DstParamData%TerrainSliceValid)) then + allocate(DstParamData%TerrainSliceValid(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%TerrainSliceValid.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%TerrainSliceValid = SrcParamData%TerrainSliceValid + end if + if (allocated(SrcParamData%TerrainSliceSeriesUn)) then + LB(1:1) = lbound(SrcParamData%TerrainSliceSeriesUn) + UB(1:1) = ubound(SrcParamData%TerrainSliceSeriesUn) + if (.not. allocated(DstParamData%TerrainSliceSeriesUn)) then + allocate(DstParamData%TerrainSliceSeriesUn(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%TerrainSliceSeriesUn.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%TerrainSliceSeriesUn = SrcParamData%TerrainSliceSeriesUn + end if + if (allocated(SrcParamData%TerrainSliceNPtsTotal)) then + LB(1:1) = lbound(SrcParamData%TerrainSliceNPtsTotal) + UB(1:1) = ubound(SrcParamData%TerrainSliceNPtsTotal) + if (.not. allocated(DstParamData%TerrainSliceNPtsTotal)) then + allocate(DstParamData%TerrainSliceNPtsTotal(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%TerrainSliceNPtsTotal.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%TerrainSliceNPtsTotal = SrcParamData%TerrainSliceNPtsTotal + end if + if (allocated(SrcParamData%TerrainSlicePtsIdx)) then + LB(1:1) = lbound(SrcParamData%TerrainSlicePtsIdx) + UB(1:1) = ubound(SrcParamData%TerrainSlicePtsIdx) + if (.not. allocated(DstParamData%TerrainSlicePtsIdx)) then + allocate(DstParamData%TerrainSlicePtsIdx(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%TerrainSlicePtsIdx.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%TerrainSlicePtsIdx = SrcParamData%TerrainSlicePtsIdx + end if + if (allocated(SrcParamData%TerrainSlicePtsFlat)) then + LB(1:2) = lbound(SrcParamData%TerrainSlicePtsFlat) + UB(1:2) = ubound(SrcParamData%TerrainSlicePtsFlat) + if (.not. allocated(DstParamData%TerrainSlicePtsFlat)) then + allocate(DstParamData%TerrainSlicePtsFlat(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%TerrainSlicePtsFlat.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%TerrainSlicePtsFlat = SrcParamData%TerrainSlicePtsFlat + 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 @@ -2182,8 +2717,8 @@ subroutine AWAE_DestroyParam(ParamData, ErrStat, ErrMsg) type(AWAE_ParameterType), intent(inout) :: ParamData integer(IntKi), intent( out) :: ErrStat character(*), intent( out) :: ErrMsg - integer(B4Ki) :: i1 - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: i1, i2 + integer(B4Ki) :: LB(2), UB(2) integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 character(*), parameter :: RoutineName = 'AWAE_DestroyParam' @@ -2233,6 +2768,63 @@ subroutine AWAE_DestroyParam(ParamData, ErrStat, ErrMsg) if (allocated(ParamData%OutDisWindYvalid)) then deallocate(ParamData%OutDisWindYvalid) end if + if (allocated(ParamData%PlaneSliceName)) then + deallocate(ParamData%PlaneSliceName) + end if + if (allocated(ParamData%PlaneSliceOrigin)) then + deallocate(ParamData%PlaneSliceOrigin) + end if + if (allocated(ParamData%PlaneSliceNormalAxis)) then + deallocate(ParamData%PlaneSliceNormalAxis) + end if + if (allocated(ParamData%PlaneSliceAxis1)) then + deallocate(ParamData%PlaneSliceAxis1) + end if + if (allocated(ParamData%PlaneSliceAxis2)) then + deallocate(ParamData%PlaneSliceAxis2) + end if + if (allocated(ParamData%PlaneSliceN1)) then + deallocate(ParamData%PlaneSliceN1) + end if + if (allocated(ParamData%PlaneSliceN2)) then + deallocate(ParamData%PlaneSliceN2) + end if + if (allocated(ParamData%PlaneSliceILoReq)) then + deallocate(ParamData%PlaneSliceILoReq) + end if + if (allocated(ParamData%PlaneSliceIHiReq)) then + deallocate(ParamData%PlaneSliceIHiReq) + end if + if (allocated(ParamData%PlaneSliceJLoReq)) then + deallocate(ParamData%PlaneSliceJLoReq) + end if + if (allocated(ParamData%PlaneSliceJHiReq)) then + deallocate(ParamData%PlaneSliceJHiReq) + end if + if (allocated(ParamData%PlaneSliceValid)) then + deallocate(ParamData%PlaneSliceValid) + end if + if (allocated(ParamData%PlaneSliceSeriesUn)) then + deallocate(ParamData%PlaneSliceSeriesUn) + end if + if (allocated(ParamData%TerrainSliceName)) then + deallocate(ParamData%TerrainSliceName) + end if + if (allocated(ParamData%TerrainSliceValid)) then + deallocate(ParamData%TerrainSliceValid) + end if + if (allocated(ParamData%TerrainSliceSeriesUn)) then + deallocate(ParamData%TerrainSliceSeriesUn) + end if + if (allocated(ParamData%TerrainSliceNPtsTotal)) then + deallocate(ParamData%TerrainSliceNPtsTotal) + end if + if (allocated(ParamData%TerrainSlicePtsIdx)) then + deallocate(ParamData%TerrainSlicePtsIdx) + end if + if (allocated(ParamData%TerrainSlicePtsFlat)) then + deallocate(ParamData%TerrainSlicePtsFlat) + end if nullify(ParamData%WAT_FlowField) end subroutine @@ -2240,8 +2832,8 @@ subroutine AWAE_PackParam(RF, Indata) type(RegFile), intent(inout) :: RF type(AWAE_ParameterType), intent(in) :: InData character(*), parameter :: RoutineName = 'AWAE_PackParam' - integer(B4Ki) :: i1 - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: i1, i2 + integer(B4Ki) :: LB(2), UB(2) logical :: PtrInIndex if (RF%ErrStat >= AbortErrLev) return call RegPack(RF, InData%WindFilePath) @@ -2297,9 +2889,36 @@ subroutine AWAE_PackParam(RF, Indata) call RegPack(RF, InData%NOutDisWindXZ) call RegPackAlloc(RF, InData%OutDisWindY) call RegPackAlloc(RF, InData%OutDisWindYvalid) + call RegPack(RF, InData%WrPlaneSkp) + call RegPack(RF, InData%NumPlaneSlices) + call RegPackAlloc(RF, InData%PlaneSliceName) + call RegPackAlloc(RF, InData%PlaneSliceOrigin) + call RegPackAlloc(RF, InData%PlaneSliceNormalAxis) + call RegPackAlloc(RF, InData%PlaneSliceAxis1) + call RegPackAlloc(RF, InData%PlaneSliceAxis2) + call RegPackAlloc(RF, InData%PlaneSliceN1) + call RegPackAlloc(RF, InData%PlaneSliceN2) + call RegPackAlloc(RF, InData%PlaneSliceILoReq) + call RegPackAlloc(RF, InData%PlaneSliceIHiReq) + call RegPackAlloc(RF, InData%PlaneSliceJLoReq) + call RegPackAlloc(RF, InData%PlaneSliceJHiReq) + call RegPackAlloc(RF, InData%PlaneSliceValid) + call RegPackAlloc(RF, InData%PlaneSliceSeriesUn) + call RegPack(RF, InData%WrTerrainSkp) + call RegPack(RF, InData%NumTerrainSlices) + call RegPackAlloc(RF, InData%TerrainSliceName) + call RegPackAlloc(RF, InData%TerrainSliceValid) + call RegPackAlloc(RF, InData%TerrainSliceSeriesUn) + call RegPackAlloc(RF, InData%TerrainSliceNPtsTotal) + call RegPackAlloc(RF, InData%TerrainSlicePtsIdx) + call RegPackAlloc(RF, InData%TerrainSlicePtsFlat) 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 @@ -2315,8 +2934,8 @@ subroutine AWAE_UnPackParam(RF, OutData) type(RegFile), intent(inout) :: RF type(AWAE_ParameterType), intent(inout) :: OutData character(*), parameter :: RoutineName = 'AWAE_UnPackParam' - integer(B4Ki) :: i1 - integer(B4Ki) :: LB(1), UB(1) + integer(B4Ki) :: i1, i2 + integer(B4Ki) :: LB(2), UB(2) integer(IntKi) :: stat logical :: IsAllocAssoc integer(B8Ki) :: PtrIdx @@ -2383,9 +3002,36 @@ subroutine AWAE_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%NOutDisWindXZ); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%OutDisWindY); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%OutDisWindYvalid); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%WrPlaneSkp); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%NumPlaneSlices); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceName); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceOrigin); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceNormalAxis); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceAxis1); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceAxis2); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceN1); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceN2); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceILoReq); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceIHiReq); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceJLoReq); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceJHiReq); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceValid); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PlaneSliceSeriesUn); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%WrTerrainSkp); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%NumTerrainSlices); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSliceName); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSliceValid); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSliceSeriesUn); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSliceNPtsTotal); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSlicePtsIdx); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TerrainSlicePtsFlat); 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..ca8df71867 --- /dev/null +++ b/modules/awae/src/AWAE_vtk.f90 @@ -0,0 +1,508 @@ +!********************************************************************************************************************************** +! 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 what plane this was first written out at (initialized to huge, so this will catch only the first) + if (m%WakeVTK_StartN(np_wp,nt_wp) > n) m%WakeVTK_StartN(np_wp,nt_wp) = n + 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 + if (m%WakeVTK_StartN(np_wp, nt_wp) > n_final) 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/tests/test_AMReX_reader.f90 b/modules/awae/tests/test_AMReX_reader.f90 index b382ed4750..37128ce7d4 100644 --- a/modules/awae/tests/test_AMReX_reader.f90 +++ b/modules/awae/tests/test_AMReX_reader.f90 @@ -27,7 +27,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 +36,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 +51,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 +110,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 +119,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 +134,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) diff --git a/modules/wakedynamics/src/WakeDynamics.f90 b/modules/wakedynamics/src/WakeDynamics.f90 index 05b357088e..b2d88e58a4 100644 --- a/modules/wakedynamics/src/WakeDynamics.f90 +++ b/modules/wakedynamics/src/WakeDynamics.f90 @@ -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 @@ -574,12 +579,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 +603,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 !............................................................................................ @@ -950,90 +966,123 @@ 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 - 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 consecutive out-of-bounds planes that are within 2*dr of each other + ! TODO: this assumes that only sequential planes will be out of bounds. We may need to modify this for drones. + ! -------------------------------------------------------------------------------- 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 + i = maxPln + do while (i >= 1) + ! Check if plane i is out of domain in any dimension + if (PlaneOutOfBounds(xd%p_plane(:,i))) then + ! Check if the adjacent lower-index plane (i-1) is also out of domain + if (PlaneOutOfBounds(xd%p_plane(:,i-1))) then + ! Check spatial proximity + if (TwoNorm(xd%p_plane(:,i) - xd%p_plane(:,i-1)) <= 2.0_ReKi * p%dr) then + call MergeWakePlanes(i-1, i) 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 - end if end if + i = i - 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 + subroutine updateVelocityPolar() integer(intKi) :: i,j real(ReKi) :: dx, absdx @@ -1123,8 +1172,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) @@ -1616,7 +1665,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 +1690,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 +1722,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 +1753,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 +1944,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..98c588a0d6 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" - @@ -186,6 +187,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..5ed41f340c 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 ======= @@ -208,6 +209,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 +459,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 +484,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 +497,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) @@ -1409,7 +1414,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 +1492,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 +1569,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 +1588,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 +1632,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/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index efe0587cab..a53fc66fa7 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -440,6 +440,8 @@ if(BUILD_FASTFARM) ff_regression("AMReX" "" "fastfarm") ff_regression("TSinflow" "" "fastfarm") ff_regression("LESinflow" "" "fastfarm") + ff_regression("LESinflow_ArbSlice" "" "fastfarm") + ff_regression("LESinflow_TerrainSlice" "" "fastfarm") ff_regression("TSinflow_curl" "" "fastfarm") ff_regression("ModAmb_3" "" "fastfarm") ff_regression("TSinflowADskSED" "" "fastfarm;aerodisk;simple-elastodyn") diff --git a/reg_tests/r-test b/reg_tests/r-test index 401b7fff63..be90f035d9 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 401b7fff63108a5b50ca824149183fe3609203a2 +Subproject commit be90f035d9a6d0de87425efdab4e062c17cfcce6 diff --git a/vs-build/modules/AWAE.vfproj b/vs-build/modules/AWAE.vfproj index 95ceae90ea..7e6227167b 100644 --- a/vs-build/modules/AWAE.vfproj +++ b/vs-build/modules/AWAE.vfproj @@ -94,6 +94,7 @@ +