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/CMakeLists.txt b/modules/aerodyn/CMakeLists.txt index 0c9fcda1a5..de74705881 100644 --- a/modules/aerodyn/CMakeLists.txt +++ b/modules/aerodyn/CMakeLists.txt @@ -74,6 +74,17 @@ 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; 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 + PROPERTIES COMPILE_FLAGS "-fno-trapping-math" + ) +endif() + # ADI lib add_library(adilib STATIC src/AeroDyn_Inflow.f90 diff --git a/modules/aerodyn/src/FVW_BiotSavart.f90 b/modules/aerodyn/src/FVW_BiotSavart.f90 index bf93782203..e73af674ae 100644 --- a/modules/aerodyn/src/FVW_BiotSavart.f90 +++ b/modules/aerodyn/src/FVW_BiotSavart.f90 @@ -364,22 +364,92 @@ 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 + ! 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 + 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 diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index b8bfb615d6..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 @@ -198,6 +198,13 @@ 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. + 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 if (Check( Inp%nNWPanels<1 , 'Number of near wake panels (`nNWPanels`) must be >=1 when using circulation solving with polar data (`CircSolvMethod=1`)')) return diff --git a/modules/aerodyn/src/FVW_VortexTools.f90 b/modules/aerodyn/src/FVW_VortexTools.f90 index 18073d8901..d30b5409cd 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 @@ -68,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 @@ -77,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 @@ -587,6 +595,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) @@ -624,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) @@ -720,6 +772,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 @@ -1279,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 @@ -1448,7 +1510,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 0) then + 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 @@ -1651,6 +1699,81 @@ recursive subroutine ui_tree_part_11(node, CP, Uind) end subroutine ui_tree_part_11 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 !< @@ -1696,6 +1819,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