From 58d6ba68a330ae2861ba18f6ac740a9e16d84da6 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Thu, 10 Sep 2026 16:39:34 -0600 Subject: [PATCH 1/8] FVW: vectorize near-field particle sum in ui_part_nograd Hoist RegFunction select outside the CP/particle loops and add masked/branchless SIMD variants for idRegNone and idRegCompact. Use r2s=max(r2,MINNORM2) floor to keep the divide finite (no Inf/NaN), with an msk gate for near-coincident points. idRegExp keeps the guarded r>2rc branch that skips the exp (exp blocks SIMD anyway). Scalar ui_part_nograd_11 kernel is unchanged. Co-authored-by: Copilot <198982749+Copilot@users.noreply.github.com> Co-authored-by: Claude --- modules/aerodyn/src/FVW_BiotSavart.f90 | 101 +++++++++++++++++++++---- 1 file changed, 85 insertions(+), 16 deletions(-) diff --git a/modules/aerodyn/src/FVW_BiotSavart.f90 b/modules/aerodyn/src/FVW_BiotSavart.f90 index bf93782203..3e67f2aec4 100644 --- a/modules/aerodyn/src/FVW_BiotSavart.f90 +++ b/modules/aerodyn/src/FVW_BiotSavart.f90 @@ -364,22 +364,91 @@ subroutine ui_part_nograd(nCPS, CPs, nPart, Part, Alpha, RegFunction, RegParam, real(ReKi), dimension(:,:), intent(in) :: Alpha !< Particle intensity [m^3/s] (3 x nPart+) omega dV= alpha integer(IntKi), intent(in) :: RegFunction !< Regularization function real(ReKi), dimension(:), intent(in) :: RegParam !< Regularization parameter (nPart+) - real(ReKi), dimension(3) :: UItmp !< - real(ReKi), dimension(3) :: DP !< - integer :: icp,ip - ! TODO: inlining of regularization - !$OMP PARALLEL DEFAULT(SHARED) - !$OMP DO PRIVATE(icp,ip, DP, UItmp) schedule(runtime) - do icp=1,nCPs ! loop on CPs - do ip=1,nPart ! loop on particles - UItmp(1:3) = 0.0_ReKi - DP(1:3) = CPs(1:3,icp)-Part(1:3,ip) - call ui_part_nograd_11(DP, Alpha(1:3,ip), RegFunction , RegParam(ip), UItmp) - UIout(1:3,icp)=UIout(1:3,icp)+UItmp(1:3) - enddo! loop on particles - enddo ! loop CPs - !$OMP END DO - !$OMP END PARALLEL + integer :: icp, ip + real(ReKi) :: CPx, CPy, CPz !< Control point coordinates (loop-invariant over particles) + real(ReKi) :: dx, dy, dz !< CP - particle + real(ReKi) :: r2, r2s, rn, r3 !< |r|^2, guarded |r|^2, |r|, |r|^3 + real(ReKi) :: msk, sp !< singularity mask (0/1) and ScalarPart + real(ReKi) :: Cx, Cy, Cz !< Alpha x r + real(ReKi) :: rc2, rc3, tt !< core^2, core^3, (r/rc)^2 + real(ReKi) :: ux, uy, uz !< per-CP accumulator (enables SIMD reduction) + real(ReKi), parameter :: MINNORM2 = MINNORM*MINNORM !< r^2 singularity guard: mask+max avoid the sqrt/branch and any 1/0 + ! Branchless, RegFunction-hoisted inner loops so the particle sum vectorizes (kernel math matches ui_part_nograd_11) + select case (RegFunction) + case (idRegNone) ! No mollification + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(icp,ip,CPx,CPy,CPz,dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz,ux,uy,uz) schedule(runtime) + do icp=1,nCPs + CPx=CPs(1,icp); CPy=CPs(2,icp); CPz=CPs(3,icp) + ux=0.0_ReKi; uy=0.0_ReKi; uz=0.0_ReKi + !$OMP SIMD reduction(+:ux,uy,uz) private(dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz) + do ip=1,nPart + dx=CPx-Part(1,ip); dy=CPy-Part(2,ip); dz=CPz-Part(3,ip) + r2 = dx*dx + dy*dy + dz*dz + msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) + r2s = max(r2, MINNORM2) + rn = sqrt(r2s); r3 = r2s*rn + Cx = Alpha(2,ip)*dz - Alpha(3,ip)*dy + Cy = Alpha(3,ip)*dx - Alpha(1,ip)*dz + Cz = Alpha(1,ip)*dy - Alpha(2,ip)*dx + sp = msk*fourpi_inv/r3 + ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp + end do + UIout(1,icp)=UIout(1,icp)+ux; UIout(2,icp)=UIout(2,icp)+uy; UIout(3,icp)=UIout(3,icp)+uz + end do + !$OMP END PARALLEL DO + case (idRegExp) ! Exp mollifier: exp() blocks SIMD anyway, so keep the r>2rc branch that skips the exp + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(icp,ip,CPx,CPy,CPz,dx,dy,dz,r2,rn,r3,sp,Cx,Cy,Cz,rc3,ux,uy,uz) schedule(runtime) + do icp=1,nCPs + CPx=CPs(1,icp); CPy=CPs(2,icp); CPz=CPs(3,icp) + ux=0.0_ReKi; uy=0.0_ReKi; uz=0.0_ReKi + do ip=1,nPart + dx=CPx-Part(1,ip); dy=CPy-Part(2,ip); dz=CPz-Part(3,ip) + r2 = dx*dx + dy*dy + dz*dz + if (r2 < MINNORM2) cycle ! on singularity + rn = sqrt(r2); r3 = r2*rn + rc3 = RegParam(ip)*RegParam(ip)*RegParam(ip) + Cx = Alpha(2,ip)*dz - Alpha(3,ip)*dy + Cy = Alpha(3,ip)*dx - Alpha(1,ip)*dz + Cz = Alpha(1,ip)*dy - Alpha(2,ip)*dx + if (r3 > PART_REG_CUT3*rc3) then ! r>2rc: mollifier->1, skip the expensive exp + sp = fourpi_inv/r3 + else + sp = (1.0_ReKi-exp(-r3/rc3))*fourpi_inv/r3 + end if + ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp + end do + UIout(1,icp)=UIout(1,icp)+ux; UIout(2,icp)=UIout(2,icp)+uy; UIout(3,icp)=UIout(3,icp)+uz + end do + !$OMP END PARALLEL DO + case (idRegCompact) ! Truly compact C2 core: exactly singular for r>=rc, rc=PART_REG_C2*RegParam + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(icp,ip,CPx,CPy,CPz,dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz,rc2,rc3,tt,ux,uy,uz) schedule(runtime) + do icp=1,nCPs + CPx=CPs(1,icp); CPy=CPs(2,icp); CPz=CPs(3,icp) + ux=0.0_ReKi; uy=0.0_ReKi; uz=0.0_ReKi + !$OMP SIMD reduction(+:ux,uy,uz) private(dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz,rc2,rc3,tt) + do ip=1,nPart + dx=CPx-Part(1,ip); dy=CPy-Part(2,ip); dz=CPz-Part(3,ip) + r2 = dx*dx + dy*dy + dz*dz + msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) + r2s = max(r2, MINNORM2) + rn = sqrt(r2s); r3 = r2s*rn + rc2 = (PART_REG_C2*RegParam(ip))**2 + rc3 = rc2*PART_REG_C2*RegParam(ip) + tt = r2s/rc2 + Cx = Alpha(2,ip)*dz - Alpha(3,ip)*dy + Cy = Alpha(3,ip)*dx - Alpha(1,ip)*dz + Cz = Alpha(1,ip)*dy - Alpha(2,ip)*dx + sp = merge(fourpi_inv/r3, (35._ReKi + tt*(-42._ReKi + 15._ReKi*tt))*fourpi_inv/(8._ReKi*rc3), r2s >= rc2) + sp = msk*sp + ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp + end do + UIout(1,icp)=UIout(1,icp)+ux; UIout(2,icp)=UIout(2,icp)+uy; UIout(3,icp)=UIout(3,icp)+uz + end do + !$OMP END PARALLEL DO + case default + print*,'[ERROR] Wrong regularization function for particles',RegFunction + STOP + end select end subroutine ui_part_nograd !> Induced velocity from 1 particle at 1 control point. The velocity gradient is not computed From 8e777f93d913d2c71749bd06fa3eec315503e5a4 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Fri, 11 Sep 2026 11:40:29 -0600 Subject: [PATCH 2/8] =?UTF-8?q?=EF=BB=BFFVW:=20give=20treecode=20leaf=20no?= =?UTF-8?q?des=20multipoles=20+=20SIMD=20leaf=20sum?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Turn childless nodes into standard treecode leaves that carry moments and are MAC-tested, so a far control point uses the leaf multipole instead of direct- summing the bucket. Move early-leaf termination after the P2M moment loop so early leaves keep valid moments, and set the single-particle root monopole so a far MAC hit there is exact. Unify the traversal: every non-empty node computes distDirect/r, near -> direct-sum its leaves and recurse branches, far -> node multipole (held short of M2M). Add contained ui_leaves_nograd with masked/ branchless !$OMP SIMD over particles for idRegNone/idRegCompact (idRegExp left scalar). Early termination triggers on a small bucket (nPart<=Ncrit) or a sub- core cell, sizing leaves for the vector loop. Co-authored-by: Copilot <198982749+Copilot@users.noreply.github.com> Co-authored-by: Claude --- modules/aerodyn/src/FVW_VortexTools.f90 | 112 +++++++++++++++++++----- 1 file changed, 92 insertions(+), 20 deletions(-) diff --git a/modules/aerodyn/src/FVW_VortexTools.f90 b/modules/aerodyn/src/FVW_VortexTools.f90 index 18073d8901..ebc85c46e8 100644 --- a/modules/aerodyn/src/FVW_VortexTools.f90 +++ b/modules/aerodyn/src/FVW_VortexTools.f90 @@ -33,6 +33,11 @@ module FVW_VortexTools ! Tree parameters integer, parameter :: IK1 = selected_int_kind(1) ! to store particle branch number (from 1 to 8) + ! Particle-tree leaf termination (implements the "minimum cell size" TODO): once a cell is at/below the core + ! scale the singular-kernel multipole is invalid there anyway, so keep its particles as a direct leaf bucket. + ! This prunes useless deep nodes and lengthens the leaf loop so it vectorizes. Ncrit caps dense sub-core clumps. + real(ReKi), parameter :: TREE_LEAF_REGFACTOR = 1.0_ReKi !< Stop subdividing when cell radius <= this * maxRegParam + integer, parameter :: TREE_LEAF_NCRIT = 32 !< and particle count <= this (conservative; sized for SIMD leaves) integer,parameter :: M0 = 1, M1_1=2, M1_2=3, M1_3=4, M2_11=5, M2_21=6, M2_22=7, M2_31=8, M2_32=9, M2_33=10 ! For moment coefficients integer,parameter :: M0_000 = 1 integer,parameter :: M1_100 = 2 @@ -587,6 +592,7 @@ subroutine grow_tree_part(Tree, nPart, PartP, PartAlpha, PartRegFunction, PartRe node%radius=0 node%center(1:3)=Part%P(1:3,1) node%Moments=0.0_ReKi + node%Moments(1:3,M0_000)=Part%Alpha(1:3,1) ! monopole, so a far MAC hit on this childless node is exact node%maxRegParam = Part%RegParam(1) nullify(node%iPart) nullify(node%branches) @@ -720,6 +726,14 @@ subroutine grow_tree_part_substep(node, Part) node%Moments(1:3,10) = node%Moments(1:3,10) + PartAlpha*DeltaP(3)*DeltaP(3) ! j=3,k=3 end do + ! --- Early leaf termination: stop if a small bucket (nPart<=Ncrit, classic Barnes-Hut leaf -> longer SIMD loop) OR sub-core (radius<=core, bounds depth for dense/coincident clumps). Moments are kept either way so far CPs still use the node multipole (standard leaf MAC). + if (node%radius <= TREE_LEAF_REGFACTOR*node%maxRegParam .or. node%nPart <= TREE_LEAF_NCRIT) then + allocate(node%leaves(1:node%nPart)) + node%leaves(1:node%nPart) = node%iPart(1:node%nPart) + if (associated(node%iPart)) deallocate(node%iPart) + return + end if + ! --- Distributing particles to the 8 octants (based on the geometric center!) allocate (PartOctant(1:node%nPart)) npart_per_octant(1:8)=0 @@ -1448,7 +1462,8 @@ end subroutine cut_tree_segment_parallel ! --- Velocity computation ! -------------------------------------------------------------------------------- subroutine ui_tree_part(Tree, icp_end, CPs, BranchFactor, DistanceDirect, Uind, ErrStat, ErrMsg) - use FVW_BiotSavart, only: ui_part_nograd_11, PartRegFloorFactor + use FVW_BiotSavart, only: PartRegFloorFactor + use FVW_BiotSavart, only: idRegNone, idRegExp, idRegCompact, PART_REG_C2, PART_REG_CUT3, MINNORM type(T_Tree), target, intent(inout) :: Tree !< integer, intent(in ) :: icp_end !< Number of CPs to use Batched direct sum over a node's leaves (particle indices), masked/branchless so it vectorizes. + !! Kernel math mirrors ui_part_nograd_11; the per-leaf gather uses the leaves() index array. + subroutine ui_leaves_nograd(CP, leaves, Uind) + real(ReKi), dimension(3), intent(in ) :: CP + integer, dimension(:), intent(in ) :: leaves + real(ReKi), dimension(3), intent(inout) :: Uind + integer :: i, ip, nL + real(ReKi) :: dx, dy, dz, r2, r2s, rn, r3, msk, sp, Cx, Cy, Cz, rc2, rc3, tt + real(ReKi) :: ux, uy, uz + real(ReKi), parameter :: MINNORM2 = MINNORM*MINNORM + nL = size(leaves) + ux=0.0_ReKi; uy=0.0_ReKi; uz=0.0_ReKi + select case (Part%RegFunction) + case (idRegNone) ! No mollification + !$OMP SIMD reduction(+:ux,uy,uz) private(ip,dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz) + do i=1,nL + ip=leaves(i) + dx=CP(1)-Part%P(1,ip); dy=CP(2)-Part%P(2,ip); dz=CP(3)-Part%P(3,ip) + r2 = dx*dx + dy*dy + dz*dz + msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) + r2s = max(r2, MINNORM2) + rn = sqrt(r2s); r3 = r2s*rn + Cx = Part%Alpha(2,ip)*dz - Part%Alpha(3,ip)*dy + Cy = Part%Alpha(3,ip)*dx - Part%Alpha(1,ip)*dz + Cz = Part%Alpha(1,ip)*dy - Part%Alpha(2,ip)*dx + sp = msk*fourpi_inv/r3 + ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp + end do + case (idRegExp) ! Exp mollifier: exp() blocks SIMD anyway, keep the r>2rc branch that skips the exp + do i=1,nL + ip=leaves(i) + dx=CP(1)-Part%P(1,ip); dy=CP(2)-Part%P(2,ip); dz=CP(3)-Part%P(3,ip) + r2 = dx*dx + dy*dy + dz*dz + if (r2 < MINNORM2) cycle ! on singularity + rn = sqrt(r2); r3 = r2*rn + rc3 = Part%RegParam(ip)*Part%RegParam(ip)*Part%RegParam(ip) + Cx = Part%Alpha(2,ip)*dz - Part%Alpha(3,ip)*dy + Cy = Part%Alpha(3,ip)*dx - Part%Alpha(1,ip)*dz + Cz = Part%Alpha(1,ip)*dy - Part%Alpha(2,ip)*dx + if (r3 > PART_REG_CUT3*rc3) then ! r>2rc: mollifier->1, skip the expensive exp + sp = fourpi_inv/r3 + else + sp = (1.0_ReKi-exp(-r3/rc3))*fourpi_inv/r3 + end if + ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp + end do + case (idRegCompact) ! Truly compact C2 core: exactly singular for r>=rc, rc=PART_REG_C2*RegParam + !$OMP SIMD reduction(+:ux,uy,uz) private(ip,dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz,rc2,rc3,tt) + do i=1,nL + ip=leaves(i) + dx=CP(1)-Part%P(1,ip); dy=CP(2)-Part%P(2,ip); dz=CP(3)-Part%P(3,ip) + r2 = dx*dx + dy*dy + dz*dz + msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) + r2s = max(r2, MINNORM2) + rn = sqrt(r2s); r3 = r2s*rn + rc2 = (PART_REG_C2*Part%RegParam(ip))**2 + rc3 = rc2*PART_REG_C2*Part%RegParam(ip) + tt = r2s/rc2 + Cx = Part%Alpha(2,ip)*dz - Part%Alpha(3,ip)*dy + Cy = Part%Alpha(3,ip)*dx - Part%Alpha(1,ip)*dz + Cz = Part%Alpha(1,ip)*dy - Part%Alpha(2,ip)*dx + sp = merge(fourpi_inv/r3, (35._ReKi + tt*(-42._ReKi + 15._ReKi*tt))*fourpi_inv/(8._ReKi*rc3), r2s >= rc2) + sp = msk*sp + ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp + end do + case default + print*,'[ERROR] Wrong regularization function for particles',Part%RegFunction + STOP + end select + Uind(1)=Uind(1)+ux; Uind(2)=Uind(2)+uy; Uind(3)=Uind(3)+uz + end subroutine ui_leaves_nograd end subroutine ui_tree_part subroutine ui_tree_segment(Tree, CPs, icp_end, BranchFactor, DistanceDirect, Uind, ErrStat, ErrMsg) From 2223d41d8430bd6fa5dd0f9473ea8aded67e402c Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Tue, 15 Sep 2026 11:55:59 -0600 Subject: [PATCH 3/8] FVW: validate TreeBranchFactor >= 1 for tree velocity methods The multipole expansion only converges when the target lies outside the source radius (BranchFactor >= 1); a value < 1 makes the series diverge. This was previously unchecked, so any value (including negative) was accepted. Enforce the bound at input time, only when a tree velocity method is selected. Co-authored-by: Copilot <198982749+Copilot@users.noreply.github.com> Co-authored-by: Claude --- modules/aerodyn/src/FVW_IO.f90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 40b03d1654..e0bd3ef232 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -197,6 +197,11 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) if (Check(.not.(ANY(idVelocityVALID ==Inp%VelocityMethod(1))), 'Velocity method (VelocityMethod(1)) not valid: '//trim(Num2LStr(Inp%VelocityMethod(1))))) return if (Check(.not.(ANY(idVelocityVALID ==Inp%VelocityMethod(2))), 'Velocity method (VelocityMethod(2)) not valid: '//trim(Num2LStr(Inp%VelocityMethod(2))))) return + ! Necessary condition for multipole convergence: BranchFactor<1 puts targets inside the source radius (series diverges). Default 1.5 adds geometric margin. + if (Inp%VelocityMethod(1)==idVelocityTreePart .or. Inp%VelocityMethod(1)==idVelocityTreeSeg) then + if (Check( Inp%TreeBranchFactor(1)<1.0_ReKi , 'Tree branch factor (TreeBranchFactor) must be >=1.')) return + endif + if (Check( Inp%DTfvw < p%DTaero, 'DTfvw must be >= DTaero from AD15.')) return if (Inp%CircSolvMethod == idCircPolarData) then if (Check( Inp%nNWPanels<1 , 'Number of near wake panels (`nNWPanels`) must be >=1 when using circulation solving with polar data (`CircSolvMethod=1`)')) return From e294c6b32645d5488f964789aafe6f8920c9109c Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Fri, 25 Sep 2026 17:31:44 -0600 Subject: [PATCH 4/8] AD OLAF: FP-hygiene floor in compact kernels, TreeBranchFactor(2) validation, segment-tree invariant note C1: floor the compact-C2 divisors (rc2/rc3) and clamp tt in the branchless idRegCompact SIMD kernels (ui_part_nograd, ui_leaves_nograd) so the always-evaluated, discarded merge arm never divides by zero or overflows when RegParam==0. Behavior-preserving. C2: validate TreeBranchFactor for both VelocityMethod indices (loop over 1:2) and document the >=1 requirement in the OLAF input-file docs. C3: document the segment-tree invariant in ui_tree_segment_11 (comment only). Co-authored-by: GitHub Copilot Co-authored-by: Claude --- docs/source/user/aerodyn-olaf/InputFiles.rst | 2 ++ modules/aerodyn/src/FVW_BiotSavart.f90 | 7 ++++--- modules/aerodyn/src/FVW_IO.f90 | 10 ++++++---- modules/aerodyn/src/FVW_VortexTools.f90 | 9 ++++++--- 4 files changed, 18 insertions(+), 10 deletions(-) diff --git a/docs/source/user/aerodyn-olaf/InputFiles.rst b/docs/source/user/aerodyn-olaf/InputFiles.rst index 19c97b188e..ee2bdff4be 100644 --- a/docs/source/user/aerodyn-olaf/InputFiles.rst +++ b/docs/source/user/aerodyn-olaf/InputFiles.rst @@ -260,6 +260,8 @@ The default option is *[2]*. **TreeBranchFactor** [-] specifies the dimensionless distance, in branch radius, above which a multipole calculation is used instead of a direct evaluation. Only used when *VelocityMethod* = *[2,4]*. +The value must be :math:`\geq 1`: a factor below 1 places target points inside the +source radius, where the multipole series diverges. Default value: 1.5. **PartPerSegment** [-] specifies the number of particles that are used when a diff --git a/modules/aerodyn/src/FVW_BiotSavart.f90 b/modules/aerodyn/src/FVW_BiotSavart.f90 index 3e67f2aec4..e73af674ae 100644 --- a/modules/aerodyn/src/FVW_BiotSavart.f90 +++ b/modules/aerodyn/src/FVW_BiotSavart.f90 @@ -432,9 +432,10 @@ subroutine ui_part_nograd(nCPS, CPs, nPart, Part, Alpha, RegFunction, RegParam, msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) r2s = max(r2, MINNORM2) rn = sqrt(r2s); r3 = r2s*rn - rc2 = (PART_REG_C2*RegParam(ip))**2 - rc3 = rc2*PART_REG_C2*RegParam(ip) - tt = r2s/rc2 + ! Floor divisors/clamp tt: merge evaluates both arms, so the discarded polynomial arm must not divide by zero when RegParam(ip)==0 + rc2 = max((PART_REG_C2*RegParam(ip))**2, MINNORM2) + rc3 = rc2*max(PART_REG_C2*RegParam(ip), MINNORM) + tt = min(r2s/rc2, 1.0_ReKi) Cx = Alpha(2,ip)*dz - Alpha(3,ip)*dy Cy = Alpha(3,ip)*dx - Alpha(1,ip)*dz Cz = Alpha(1,ip)*dy - Alpha(2,ip)*dx diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 025eb9ee8f..b540c7067b 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -19,7 +19,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) ! Local variables character(1024) :: PriPath ! the path to the primary input file character(1024) :: sDummy, sLine ! string to temporarially hold value of read line - integer(IntKi) :: UnIn, i + integer(IntKi) :: UnIn, i, iVel integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 ErrStat = ErrID_None @@ -199,9 +199,11 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) if (Check(.not.(ANY(idVelocityVALID ==Inp%VelocityMethod(2))), 'Velocity method (VelocityMethod(2)) not valid: '//trim(Num2LStr(Inp%VelocityMethod(2))))) return ! Necessary condition for multipole convergence: BranchFactor<1 puts targets inside the source radius (series diverges). Default 1.5 adds geometric margin. - if (Inp%VelocityMethod(1)==idVelocityTreePart .or. Inp%VelocityMethod(1)==idVelocityTreeSeg) then - if (Check( Inp%TreeBranchFactor(1)<1.0_ReKi , 'Tree branch factor (TreeBranchFactor) must be >=1.')) return - endif + do iVel = 1,2 + if (Inp%VelocityMethod(iVel)==idVelocityTreePart .or. Inp%VelocityMethod(iVel)==idVelocityTreeSeg) then + if (Check( Inp%TreeBranchFactor(iVel)<1.0_ReKi , 'Tree branch factor (TreeBranchFactor('//trim(Num2LStr(iVel))//')) must be >=1.')) return + endif + enddo if (Check( Inp%DTfvw < p%DTaero, 'DTfvw must be >= DTaero from AD15.')) return if (Inp%CircSolvMethod == idCircPolarData) then diff --git a/modules/aerodyn/src/FVW_VortexTools.f90 b/modules/aerodyn/src/FVW_VortexTools.f90 index ebc85c46e8..f2afd172c0 100644 --- a/modules/aerodyn/src/FVW_VortexTools.f90 +++ b/modules/aerodyn/src/FVW_VortexTools.f90 @@ -1705,9 +1705,10 @@ subroutine ui_leaves_nograd(CP, leaves, Uind) msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) r2s = max(r2, MINNORM2) rn = sqrt(r2s); r3 = r2s*rn - rc2 = (PART_REG_C2*Part%RegParam(ip))**2 - rc3 = rc2*PART_REG_C2*Part%RegParam(ip) - tt = r2s/rc2 + ! Floor divisors/clamp tt: merge evaluates both arms, so the discarded polynomial arm must not divide by zero when RegParam(ip)==0 + rc2 = max((PART_REG_C2*Part%RegParam(ip))**2, MINNORM2) + rc3 = rc2*max(PART_REG_C2*Part%RegParam(ip), MINNORM) + tt = min(r2s/rc2, 1.0_ReKi) Cx = Part%Alpha(2,ip)*dz - Part%Alpha(3,ip)*dy Cy = Part%Alpha(3,ip)*dx - Part%Alpha(1,ip)*dz Cz = Part%Alpha(1,ip)*dy - Part%Alpha(2,ip)*dx @@ -1768,6 +1769,8 @@ recursive subroutine ui_tree_segment_11(node, CP, Uind) integer :: iPart if (node%nPart<=0) then ! We skip the dead leaf + ! Invariant: a childless node here never carries valid moments, because segment-tree growth has no leaf-termination block (unlike ui_tree_part_11). + ! If leaf termination is ever added to grow_tree_segment_*, this must switch to the MAC-on-every-node structure of ui_tree_part_11, else far control points silently get zero velocity. elseif (.not.associated(node%branches)) then ! Loop on leaves if(associated(node%leaves)) then From dd4f2c9607db36ce2b66ec80e26775281a9073bd Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Fri, 25 Sep 2026 18:07:34 -0600 Subject: [PATCH 5/8] AD/OLAF: scope -fno-trapping-math to FVW kernels to unlock SIMD gfortran's default -ftrapping-math blocks if-conversion of the branchless compact-C2 merge in the OLAF vortex kernels, preventing vectorization. The C1 divisor floors make disabling trapping-math safe (discarded merge arm can't divide by zero). Scoped to FVW_BiotSavart.f90 and FVW_VortexTools.f90, and skipped when FPE_TRAP_ENABLED so trap-checking builds keep meaningful FPE traps. Co-authored-by: GitHub Copilot Co-authored-by: Claude --- modules/aerodyn/CMakeLists.txt | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/modules/aerodyn/CMakeLists.txt b/modules/aerodyn/CMakeLists.txt index 0c9fcda1a5..d8fff24721 100644 --- a/modules/aerodyn/CMakeLists.txt +++ b/modules/aerodyn/CMakeLists.txt @@ -74,6 +74,16 @@ add_library(aerodynlib STATIC ) target_link_libraries(aerodynlib basicaerolib aeroacousticslib seastlib nwtclibs) +# OLAF vortex kernels: gfortran's default -ftrapping-math blocks SIMD if-conversion of the branchless +# compact-C2 merge; C1's divisor floors make disabling it safe here. Skip under FPE trapping to keep that build's traps meaningful. +if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU" AND NOT FPE_TRAP_ENABLED) + set_source_files_properties( + src/FVW_BiotSavart.f90 + src/FVW_VortexTools.f90 + PROPERTIES COMPILE_FLAGS "-fno-trapping-math" + ) +endif() + # ADI lib add_library(adilib STATIC src/AeroDyn_Inflow.f90 From 404e2643bb991292a056f4d0105bb1469927ea9e Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Fri, 25 Sep 2026 18:42:35 -0600 Subject: [PATCH 6/8] AD/OLAF: pack treecode leaf sources into contiguous LeafBuf for SIMD P2P The treecode leaf kernel (ui_leaves_nograd) read sources via ip=leaves(i), an indirect stride-3 gather that blocked SIMD (idRegNone leaf failed with 'no vectype for stmt' on the SSE2 baseline). A post-build O(n) pass (compact_tree_part) copies each node's leaf particles into a contiguous, column-major Tree%LeafBuf(nAct,7) [Px Py Pz Ax Ay Az Rc], and each node records its bucket slice via iLeafStart/iLeafN. The kernel now reads unit-stride columns -> packed loads. Result-neutral (scratch copy; the shared Part arrays and the control points are untouched). Empirically both vectorizable leaf loops now vectorize (16-byte): idRegNone (previously gather-blocked) and idRegCompact; idRegExp stays scalar by design (exp control flow). All OLAF regression tests pass within tolerance. Co-authored-by: GitHub Copilot Co-authored-by: Claude --- modules/aerodyn/src/FVW_VortexTools.f90 | 110 +++++++++++++++++------- 1 file changed, 78 insertions(+), 32 deletions(-) diff --git a/modules/aerodyn/src/FVW_VortexTools.f90 b/modules/aerodyn/src/FVW_VortexTools.f90 index f2afd172c0..0b75017efd 100644 --- a/modules/aerodyn/src/FVW_VortexTools.f90 +++ b/modules/aerodyn/src/FVW_VortexTools.f90 @@ -73,6 +73,8 @@ module FVW_VortexTools integer,dimension(:),pointer :: leaves=>null() ! NOTE: leaves are introduced to save memory type(T_Node),dimension(:), pointer :: branches =>null() integer :: nPart = -1 ! Number of particles in branches and leaves of this node + integer :: iLeafStart = -1 !< First row of Tree%LeafBuf for this node's leaf bucket (-1 = no leaves) + integer :: iLeafN = 0 !< Leaf-particle count in this node (NOT nPart: a split node has leaves+branches) end type T_Node !> The type tree contains some basic data, a chained-list of nodes, and a pointer to the Particle data that were used @@ -82,6 +84,7 @@ module FVW_VortexTools integer :: iStep =-1 !< Time step at which the tree was built logical :: bGrown =.false. !< Is the tree build type(T_Node) :: Root !< Contains the chained-list of nodes + real(ReKi), dimension(:,:), allocatable :: LeafBuf !< Leaf sources packed by node bucket, cols [Px Py Pz Ax Ay Az Rc]; column-major => contiguous SIMD loads in leaf P2P (kills the leaves(i) gather) end type T_Tree interface cut_tree @@ -630,10 +633,53 @@ subroutine grow_tree_part(Tree, nPart, PartP, PartAlpha, PartRegFunction, PartRe call grow_tree_part_parallel(Tree%root, Tree%Part) ! call grow_tree_rec(Tree%root, Tree%Part) endif + ! --- Pack leaf sources contiguously for SIMD leaf P2P (kills the leaves(i) gather in ui_leaves_nograd) + if (Part%n>0) call compact_tree_part(Tree) Tree%iStep = iStep Tree%bGrown = .true. end subroutine grow_tree_part + !> Copy each node's leaf-particle sources into the contiguous Tree%LeafBuf, recording per-node [iLeafStart,iLeafN]. + !! Post-build, O(n), once per rebuild. Reorders leaf inputs only (a scratch copy); the shared Part arrays are untouched. + subroutine compact_tree_part(Tree) + type(T_Tree), intent(inout), target :: Tree + type(T_VPart), pointer :: Part + integer :: k + Part => Tree%Part + if (allocated(Tree%LeafBuf)) then + if (size(Tree%LeafBuf,1) /= Part%n) deallocate(Tree%LeafBuf) + endif + if (.not.allocated(Tree%LeafBuf)) allocate(Tree%LeafBuf(Part%n,7)) + k = 0 + call compact_node(Tree%Root) + contains + recursive subroutine compact_node(node) + type(T_Node), intent(inout) :: node + integer :: i, ip, nL + if (node%nPart<=0) return + if (associated(node%leaves)) then + nL = size(node%leaves) + node%iLeafStart = k+1 + node%iLeafN = nL + do i=1,nL + ip = node%leaves(i) + Tree%LeafBuf(k+i,1) = Part%P(1,ip); Tree%LeafBuf(k+i,2) = Part%P(2,ip); Tree%LeafBuf(k+i,3) = Part%P(3,ip) + Tree%LeafBuf(k+i,4) = Part%Alpha(1,ip); Tree%LeafBuf(k+i,5) = Part%Alpha(2,ip); Tree%LeafBuf(k+i,6) = Part%Alpha(3,ip) + Tree%LeafBuf(k+i,7) = Part%RegParam(ip) + end do + k = k + nL + else + node%iLeafStart = -1 + node%iLeafN = 0 + endif + if (associated(node%branches)) then + do i=1,size(node%branches) + call compact_node(node%branches(i)) + end do + endif + end subroutine compact_node + end subroutine compact_tree_part + !> Recursive function to grow/setup a tree. !! Note, needed preliminary calc are done by grow_tree before recursive subroutine grow_tree_part_rec(node, Part) @@ -1293,6 +1339,8 @@ subroutine cut_tree_parallel(Tree, deallocPart, deallocSgmt) logical, optional, intent(in) :: deallocPart logical, optional, intent(in) :: deallocSgmt integer :: i,i1,i2,nBranches,istat + ! --- Freeing the packed leaf buffer + if (allocated(Tree%LeafBuf)) deallocate(Tree%LeafBuf) ! --- Deallocating data we are pointing to, only if user requests it if (present(deallocPart)) then if (deallocPart) then @@ -1516,8 +1564,8 @@ recursive subroutine ui_tree_part_11(node, CP, Uind) ! Test if the control point is too close so that a direct evaluation is needed if (r0) then + call ui_leaves_nograd(CP, node%iLeafStart, node%iLeafN, Uind) endif if(associated(node%branches)) then ! TODO: consider implementing a recursive method for that: direct call on all children @@ -1650,45 +1698,44 @@ recursive subroutine ui_tree_part_11(node, CP, Uind) end if ! had more than 1 particles end subroutine ui_tree_part_11 - !> Batched direct sum over a node's leaves (particle indices), masked/branchless so it vectorizes. - !! Kernel math mirrors ui_part_nograd_11; the per-leaf gather uses the leaves() index array. - subroutine ui_leaves_nograd(CP, leaves, Uind) + !> Batched direct sum over a node's leaves, masked/branchless so it vectorizes. + !! Kernel math mirrors ui_part_nograd_11; sources read from contiguous Tree%LeafBuf rows -> packed loads (was an indirect leaves(i) stride-3 gather). + subroutine ui_leaves_nograd(CP, iStart, nL, Uind) real(ReKi), dimension(3), intent(in ) :: CP - integer, dimension(:), intent(in ) :: leaves + integer, intent(in ) :: iStart !< First LeafBuf row of this node's bucket + integer, intent(in ) :: nL !< Number of leaf particles in this node real(ReKi), dimension(3), intent(inout) :: Uind - integer :: i, ip, nL + integer :: i, iEnd real(ReKi) :: dx, dy, dz, r2, r2s, rn, r3, msk, sp, Cx, Cy, Cz, rc2, rc3, tt real(ReKi) :: ux, uy, uz real(ReKi), parameter :: MINNORM2 = MINNORM*MINNORM - nL = size(leaves) + iEnd = iStart+nL-1 ux=0.0_ReKi; uy=0.0_ReKi; uz=0.0_ReKi select case (Part%RegFunction) case (idRegNone) ! No mollification - !$OMP SIMD reduction(+:ux,uy,uz) private(ip,dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz) - do i=1,nL - ip=leaves(i) - dx=CP(1)-Part%P(1,ip); dy=CP(2)-Part%P(2,ip); dz=CP(3)-Part%P(3,ip) + !$OMP SIMD reduction(+:ux,uy,uz) private(dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz) + do i=iStart,iEnd + dx=CP(1)-Tree%LeafBuf(i,1); dy=CP(2)-Tree%LeafBuf(i,2); dz=CP(3)-Tree%LeafBuf(i,3) r2 = dx*dx + dy*dy + dz*dz msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) r2s = max(r2, MINNORM2) rn = sqrt(r2s); r3 = r2s*rn - Cx = Part%Alpha(2,ip)*dz - Part%Alpha(3,ip)*dy - Cy = Part%Alpha(3,ip)*dx - Part%Alpha(1,ip)*dz - Cz = Part%Alpha(1,ip)*dy - Part%Alpha(2,ip)*dx + Cx = Tree%LeafBuf(i,5)*dz - Tree%LeafBuf(i,6)*dy + Cy = Tree%LeafBuf(i,6)*dx - Tree%LeafBuf(i,4)*dz + Cz = Tree%LeafBuf(i,4)*dy - Tree%LeafBuf(i,5)*dx sp = msk*fourpi_inv/r3 ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp end do case (idRegExp) ! Exp mollifier: exp() blocks SIMD anyway, keep the r>2rc branch that skips the exp - do i=1,nL - ip=leaves(i) - dx=CP(1)-Part%P(1,ip); dy=CP(2)-Part%P(2,ip); dz=CP(3)-Part%P(3,ip) + do i=iStart,iEnd + dx=CP(1)-Tree%LeafBuf(i,1); dy=CP(2)-Tree%LeafBuf(i,2); dz=CP(3)-Tree%LeafBuf(i,3) r2 = dx*dx + dy*dy + dz*dz if (r2 < MINNORM2) cycle ! on singularity rn = sqrt(r2); r3 = r2*rn - rc3 = Part%RegParam(ip)*Part%RegParam(ip)*Part%RegParam(ip) - Cx = Part%Alpha(2,ip)*dz - Part%Alpha(3,ip)*dy - Cy = Part%Alpha(3,ip)*dx - Part%Alpha(1,ip)*dz - Cz = Part%Alpha(1,ip)*dy - Part%Alpha(2,ip)*dx + rc3 = Tree%LeafBuf(i,7)*Tree%LeafBuf(i,7)*Tree%LeafBuf(i,7) + Cx = Tree%LeafBuf(i,5)*dz - Tree%LeafBuf(i,6)*dy + Cy = Tree%LeafBuf(i,6)*dx - Tree%LeafBuf(i,4)*dz + Cz = Tree%LeafBuf(i,4)*dy - Tree%LeafBuf(i,5)*dx if (r3 > PART_REG_CUT3*rc3) then ! r>2rc: mollifier->1, skip the expensive exp sp = fourpi_inv/r3 else @@ -1697,21 +1744,20 @@ subroutine ui_leaves_nograd(CP, leaves, Uind) ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp end do case (idRegCompact) ! Truly compact C2 core: exactly singular for r>=rc, rc=PART_REG_C2*RegParam - !$OMP SIMD reduction(+:ux,uy,uz) private(ip,dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz,rc2,rc3,tt) - do i=1,nL - ip=leaves(i) - dx=CP(1)-Part%P(1,ip); dy=CP(2)-Part%P(2,ip); dz=CP(3)-Part%P(3,ip) + !$OMP SIMD reduction(+:ux,uy,uz) private(dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz,rc2,rc3,tt) + do i=iStart,iEnd + dx=CP(1)-Tree%LeafBuf(i,1); dy=CP(2)-Tree%LeafBuf(i,2); dz=CP(3)-Tree%LeafBuf(i,3) r2 = dx*dx + dy*dy + dz*dz msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) r2s = max(r2, MINNORM2) rn = sqrt(r2s); r3 = r2s*rn - ! Floor divisors/clamp tt: merge evaluates both arms, so the discarded polynomial arm must not divide by zero when RegParam(ip)==0 - rc2 = max((PART_REG_C2*Part%RegParam(ip))**2, MINNORM2) - rc3 = rc2*max(PART_REG_C2*Part%RegParam(ip), MINNORM) + ! Floor divisors/clamp tt: merge evaluates both arms, so the discarded polynomial arm must not divide by zero when RegParam==0 + rc2 = max((PART_REG_C2*Tree%LeafBuf(i,7))**2, MINNORM2) + rc3 = rc2*max(PART_REG_C2*Tree%LeafBuf(i,7), MINNORM) tt = min(r2s/rc2, 1.0_ReKi) - Cx = Part%Alpha(2,ip)*dz - Part%Alpha(3,ip)*dy - Cy = Part%Alpha(3,ip)*dx - Part%Alpha(1,ip)*dz - Cz = Part%Alpha(1,ip)*dy - Part%Alpha(2,ip)*dx + Cx = Tree%LeafBuf(i,5)*dz - Tree%LeafBuf(i,6)*dy + Cy = Tree%LeafBuf(i,6)*dx - Tree%LeafBuf(i,4)*dz + Cz = Tree%LeafBuf(i,4)*dy - Tree%LeafBuf(i,5)*dx sp = merge(fourpi_inv/r3, (35._ReKi + tt*(-42._ReKi + 15._ReKi*tt))*fourpi_inv/(8._ReKi*rc3), r2s >= rc2) sp = msk*sp ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp From 3d4114f53b013c0f56bfdbcc720ee36d94f8c2c0 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Fri, 25 Sep 2026 21:51:31 -0600 Subject: [PATCH 7/8] Also skip -fno-trapping-math when -ffpe-trap is set in CMAKE_Fortran_FLAGS FPE_TRAP_ENABLED only toggles the -DFPE_TRAP_ENABLED define, not any -ffpe-trap compile flag, so a user enabling traps via -DCMAKE_Fortran_FLAGS=-ffpe-trap=... would still get -fno-trapping-math on the FVW kernels. Guard on that flag too. Co-authored-by: GitHub Copilot Co-authored-by: Claude --- modules/aerodyn/CMakeLists.txt | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/modules/aerodyn/CMakeLists.txt b/modules/aerodyn/CMakeLists.txt index d8fff24721..de74705881 100644 --- a/modules/aerodyn/CMakeLists.txt +++ b/modules/aerodyn/CMakeLists.txt @@ -75,8 +75,9 @@ add_library(aerodynlib STATIC target_link_libraries(aerodynlib basicaerolib aeroacousticslib seastlib nwtclibs) # OLAF vortex kernels: gfortran's default -ftrapping-math blocks SIMD if-conversion of the branchless -# compact-C2 merge; C1's divisor floors make disabling it safe here. Skip under FPE trapping to keep that build's traps meaningful. -if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU" AND NOT FPE_TRAP_ENABLED) +# compact-C2 merge; Skip when FPE trapping is requested either via FPE_TRAP_ENABLED or an -ffpe-trap +# in CMAKE_Fortran_FLAGS, so those traps stay meaningful. +if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU" AND NOT FPE_TRAP_ENABLED AND NOT CMAKE_Fortran_FLAGS MATCHES "-ffpe-trap") set_source_files_properties( src/FVW_BiotSavart.f90 src/FVW_VortexTools.f90 From 0ca7715ab385c3ee17398d6447a94bde419fc864 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Fri, 25 Sep 2026 22:50:43 -0600 Subject: [PATCH 8/8] Promote ui_leaves_nograd to module level with contiguous LeafBuf dummy As a contained procedure of ui_tree_part, the leaf kernel host-associated Tree%LeafBuf; under -fopenmp the enclosing PARALLEL DO outlining lost the descriptor's unit-stride proof and degraded the packed loads to a gather, so the None/Compact leaf loops did not vectorize in an OpenMP build (the default VelocityMethod path). Passing LeafBuf as an explicit contiguous dummy re-establishes unit stride in the callee independent of caller outlining; both leaf loops now vectorize 16B under OPENMP=ON, with no change under OPENMP=OFF. Co-authored-by: GitHub Copilot Co-authored-by: Claude --- modules/aerodyn/src/FVW_VortexTools.f90 | 148 ++++++++++++------------ 1 file changed, 76 insertions(+), 72 deletions(-) diff --git a/modules/aerodyn/src/FVW_VortexTools.f90 b/modules/aerodyn/src/FVW_VortexTools.f90 index 0b75017efd..d30b5409cd 100644 --- a/modules/aerodyn/src/FVW_VortexTools.f90 +++ b/modules/aerodyn/src/FVW_VortexTools.f90 @@ -1565,7 +1565,7 @@ recursive subroutine ui_tree_part_11(node, CP, Uind) if (r0) then - call ui_leaves_nograd(CP, node%iLeafStart, node%iLeafN, Uind) + call ui_leaves_nograd(CP, Tree%LeafBuf, node%iLeafStart, node%iLeafN, Part%RegFunction, Uind) endif if(associated(node%branches)) then ! TODO: consider implementing a recursive method for that: direct call on all children @@ -1697,79 +1697,83 @@ recursive subroutine ui_tree_part_11(node, CP, Uind) end if ! Far enough end if ! had more than 1 particles end subroutine ui_tree_part_11 - - !> Batched direct sum over a node's leaves, masked/branchless so it vectorizes. - !! Kernel math mirrors ui_part_nograd_11; sources read from contiguous Tree%LeafBuf rows -> packed loads (was an indirect leaves(i) stride-3 gather). - subroutine ui_leaves_nograd(CP, iStart, nL, Uind) - real(ReKi), dimension(3), intent(in ) :: CP - integer, intent(in ) :: iStart !< First LeafBuf row of this node's bucket - integer, intent(in ) :: nL !< Number of leaf particles in this node - real(ReKi), dimension(3), intent(inout) :: Uind - integer :: i, iEnd - real(ReKi) :: dx, dy, dz, r2, r2s, rn, r3, msk, sp, Cx, Cy, Cz, rc2, rc3, tt - real(ReKi) :: ux, uy, uz - real(ReKi), parameter :: MINNORM2 = MINNORM*MINNORM - iEnd = iStart+nL-1 - ux=0.0_ReKi; uy=0.0_ReKi; uz=0.0_ReKi - select case (Part%RegFunction) - case (idRegNone) ! No mollification - !$OMP SIMD reduction(+:ux,uy,uz) private(dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz) - do i=iStart,iEnd - dx=CP(1)-Tree%LeafBuf(i,1); dy=CP(2)-Tree%LeafBuf(i,2); dz=CP(3)-Tree%LeafBuf(i,3) - r2 = dx*dx + dy*dy + dz*dz - msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) - r2s = max(r2, MINNORM2) - rn = sqrt(r2s); r3 = r2s*rn - Cx = Tree%LeafBuf(i,5)*dz - Tree%LeafBuf(i,6)*dy - Cy = Tree%LeafBuf(i,6)*dx - Tree%LeafBuf(i,4)*dz - Cz = Tree%LeafBuf(i,4)*dy - Tree%LeafBuf(i,5)*dx - sp = msk*fourpi_inv/r3 - ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp - end do - case (idRegExp) ! Exp mollifier: exp() blocks SIMD anyway, keep the r>2rc branch that skips the exp - do i=iStart,iEnd - dx=CP(1)-Tree%LeafBuf(i,1); dy=CP(2)-Tree%LeafBuf(i,2); dz=CP(3)-Tree%LeafBuf(i,3) - r2 = dx*dx + dy*dy + dz*dz - if (r2 < MINNORM2) cycle ! on singularity - rn = sqrt(r2); r3 = r2*rn - rc3 = Tree%LeafBuf(i,7)*Tree%LeafBuf(i,7)*Tree%LeafBuf(i,7) - Cx = Tree%LeafBuf(i,5)*dz - Tree%LeafBuf(i,6)*dy - Cy = Tree%LeafBuf(i,6)*dx - Tree%LeafBuf(i,4)*dz - Cz = Tree%LeafBuf(i,4)*dy - Tree%LeafBuf(i,5)*dx - if (r3 > PART_REG_CUT3*rc3) then ! r>2rc: mollifier->1, skip the expensive exp - sp = fourpi_inv/r3 - else - sp = (1.0_ReKi-exp(-r3/rc3))*fourpi_inv/r3 - end if - ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp - end do - case (idRegCompact) ! Truly compact C2 core: exactly singular for r>=rc, rc=PART_REG_C2*RegParam - !$OMP SIMD reduction(+:ux,uy,uz) private(dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz,rc2,rc3,tt) - do i=iStart,iEnd - dx=CP(1)-Tree%LeafBuf(i,1); dy=CP(2)-Tree%LeafBuf(i,2); dz=CP(3)-Tree%LeafBuf(i,3) - r2 = dx*dx + dy*dy + dz*dz - msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) - r2s = max(r2, MINNORM2) - rn = sqrt(r2s); r3 = r2s*rn - ! Floor divisors/clamp tt: merge evaluates both arms, so the discarded polynomial arm must not divide by zero when RegParam==0 - rc2 = max((PART_REG_C2*Tree%LeafBuf(i,7))**2, MINNORM2) - rc3 = rc2*max(PART_REG_C2*Tree%LeafBuf(i,7), MINNORM) - tt = min(r2s/rc2, 1.0_ReKi) - Cx = Tree%LeafBuf(i,5)*dz - Tree%LeafBuf(i,6)*dy - Cy = Tree%LeafBuf(i,6)*dx - Tree%LeafBuf(i,4)*dz - Cz = Tree%LeafBuf(i,4)*dy - Tree%LeafBuf(i,5)*dx - sp = merge(fourpi_inv/r3, (35._ReKi + tt*(-42._ReKi + 15._ReKi*tt))*fourpi_inv/(8._ReKi*rc3), r2s >= rc2) - sp = msk*sp - ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp - end do - case default - print*,'[ERROR] Wrong regularization function for particles',Part%RegFunction - STOP - end select - Uind(1)=Uind(1)+ux; Uind(2)=Uind(2)+uy; Uind(3)=Uind(3)+uz - end subroutine ui_leaves_nograd end subroutine ui_tree_part + !> Batched direct sum over a node's leaves, masked/branchless so it vectorizes. + !! Module-level with a contiguous LeafBuf dummy so unit stride survives OpenMP outlining of the caller; + !! as a contained procedure the host-associated buffer degraded to a gather under -fopenmp. + subroutine ui_leaves_nograd(CP, LeafBuf, iStart, nL, RegFunction, Uind) + use FVW_BiotSavart, only: idRegNone, idRegExp, idRegCompact, PART_REG_C2, PART_REG_CUT3, MINNORM + real(ReKi), dimension(3), intent(in ) :: CP + real(ReKi), dimension(:,:), contiguous, intent(in ) :: LeafBuf !< Packed leaf sources, cols [Px Py Pz Ax Ay Az Rc] + integer, intent(in ) :: iStart !< First LeafBuf row of this node's bucket + integer, intent(in ) :: nL !< Number of leaf particles in this node + integer, intent(in ) :: RegFunction + real(ReKi), dimension(3), intent(inout) :: Uind + integer :: i, iEnd + real(ReKi) :: dx, dy, dz, r2, r2s, rn, r3, msk, sp, Cx, Cy, Cz, rc2, rc3, tt + real(ReKi) :: ux, uy, uz + real(ReKi), parameter :: MINNORM2 = MINNORM*MINNORM + iEnd = iStart+nL-1 + ux=0.0_ReKi; uy=0.0_ReKi; uz=0.0_ReKi + select case (RegFunction) + case (idRegNone) ! No mollification + !$OMP SIMD reduction(+:ux,uy,uz) private(dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz) + do i=iStart,iEnd + dx=CP(1)-LeafBuf(i,1); dy=CP(2)-LeafBuf(i,2); dz=CP(3)-LeafBuf(i,3) + r2 = dx*dx + dy*dy + dz*dz + msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) + r2s = max(r2, MINNORM2) + rn = sqrt(r2s); r3 = r2s*rn + Cx = LeafBuf(i,5)*dz - LeafBuf(i,6)*dy + Cy = LeafBuf(i,6)*dx - LeafBuf(i,4)*dz + Cz = LeafBuf(i,4)*dy - LeafBuf(i,5)*dx + sp = msk*fourpi_inv/r3 + ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp + end do + case (idRegExp) ! Exp mollifier: exp() blocks SIMD anyway, keep the r>2rc branch that skips the exp + do i=iStart,iEnd + dx=CP(1)-LeafBuf(i,1); dy=CP(2)-LeafBuf(i,2); dz=CP(3)-LeafBuf(i,3) + r2 = dx*dx + dy*dy + dz*dz + if (r2 < MINNORM2) cycle ! on singularity + rn = sqrt(r2); r3 = r2*rn + rc3 = LeafBuf(i,7)*LeafBuf(i,7)*LeafBuf(i,7) + Cx = LeafBuf(i,5)*dz - LeafBuf(i,6)*dy + Cy = LeafBuf(i,6)*dx - LeafBuf(i,4)*dz + Cz = LeafBuf(i,4)*dy - LeafBuf(i,5)*dx + if (r3 > PART_REG_CUT3*rc3) then ! r>2rc: mollifier->1, skip the expensive exp + sp = fourpi_inv/r3 + else + sp = (1.0_ReKi-exp(-r3/rc3))*fourpi_inv/r3 + end if + ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp + end do + case (idRegCompact) ! Truly compact C2 core: exactly singular for r>=rc, rc=PART_REG_C2*RegParam + !$OMP SIMD reduction(+:ux,uy,uz) private(dx,dy,dz,r2,r2s,rn,r3,msk,sp,Cx,Cy,Cz,rc2,rc3,tt) + do i=iStart,iEnd + dx=CP(1)-LeafBuf(i,1); dy=CP(2)-LeafBuf(i,2); dz=CP(3)-LeafBuf(i,3) + r2 = dx*dx + dy*dy + dz*dz + msk = merge(1.0_ReKi, 0.0_ReKi, r2>=MINNORM2) + r2s = max(r2, MINNORM2) + rn = sqrt(r2s); r3 = r2s*rn + ! Floor divisors/clamp tt: merge evaluates both arms, so the discarded polynomial arm must not divide by zero when RegParam==0 + rc2 = max((PART_REG_C2*LeafBuf(i,7))**2, MINNORM2) + rc3 = rc2*max(PART_REG_C2*LeafBuf(i,7), MINNORM) + tt = min(r2s/rc2, 1.0_ReKi) + Cx = LeafBuf(i,5)*dz - LeafBuf(i,6)*dy + Cy = LeafBuf(i,6)*dx - LeafBuf(i,4)*dz + Cz = LeafBuf(i,4)*dy - LeafBuf(i,5)*dx + sp = merge(fourpi_inv/r3, (35._ReKi + tt*(-42._ReKi + 15._ReKi*tt))*fourpi_inv/(8._ReKi*rc3), r2s >= rc2) + sp = msk*sp + ux=ux+Cx*sp; uy=uy+Cy*sp; uz=uz+Cz*sp + end do + case default + print*,'[ERROR] Wrong regularization function for particles',RegFunction + STOP + end select + Uind(1)=Uind(1)+ux; Uind(2)=Uind(2)+uy; Uind(3)=Uind(3)+uz + end subroutine ui_leaves_nograd + subroutine ui_tree_segment(Tree, CPs, icp_end, BranchFactor, DistanceDirect, Uind, ErrStat, ErrMsg) use FVW_BiotSavart, only: ui_seg_11 type(T_Tree), target, intent(inout) :: Tree !<