Merge pull request #3430 from luwang00/b/OLAF

AD: fix OLAF treecode near-core regularization floor and TwrInfl OpenMP data race

Author: andrew-platt

modules/aerodyn/src/AeroDyn.f90

@@ -5605,8 +5605,9 @@ SUBROUTINE TwrInflArray( p, u, RotInflow, m, Positions, Inflow, ErrStat, ErrMsg
    ! these models are valid for only small tower deflections; check for potential division-by-zero errors:   
    call CheckTwrInfl( u, ErrStat2, ErrMsg2 ); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ); if (ErrStat >= AbortErrLev) return
 
+   ! FirstWarn_TowerStrike is firstprivate so each thread starts from the .false. set above (avoids reading uninitialized memory); ErrStat2/ErrMsg2 are private to avoid a data race
    !$OMP PARALLEL default(shared)
-   !$OMP do private(i,Pos,theta_tower_trans,W_tower,xbar,ybar,zbar,TwrCd,TwrTI,TwrClrnc,FirstWarn_TowerStrike,DisturbInflow,v) schedule(runtime)
+   !$OMP do private(i,Pos,theta_tower_trans,W_tower,xbar,ybar,zbar,TwrCd,TwrTI,TwrClrnc,DisturbInflow,v,ErrStat2,ErrMsg2) firstprivate(FirstWarn_TowerStrike) schedule(runtime)
    do i = 1, size(Positions,2)
       Pos=Positions(1:3,i)
 

modules/aerodyn/src/FVW_BiotSavart.f90

@@ -10,6 +10,8 @@ module FVW_BiotSavart
    real(ReKi),parameter :: PRECISION_UI  = epsilon(1.0_ReKi)/100 !< NOTE assuming problem of size 1
    real(ReKi),parameter :: PRECISION_EPS =  epsilon(1.0_ReKi) !< Machine Precision For the given ReKi for problems of scale 1!
    real(ReKi),parameter :: MIN_EXP_VALUE=-10.0_ReKi
+   real(ReKi),parameter :: PART_REG_NRAD =  2.0_ReKi             !< Particle exp mollifier treated as 1 beyond this many core radii (matches 2*rc far-field multipole floor)
+   real(ReKi),parameter :: PART_REG_CUT3 =  PART_REG_NRAD**3     !< Corresponding (r/rc)^3 cutoff
    real(ReKi),parameter :: MINDENOM=0.0_ReKi
 !    real(ReKi),parameter :: MINDENOM=1e-15_ReKi
    real(ReKi),parameter :: MINNORM=1e-4
@@ -374,26 +376,35 @@ subroutine ui_part_nograd_11(DeltaP, Alpha, RegFunction, RegParam, Ui)
    real(ReKi),dimension(3) :: C          !< Cross product of Alpha and r
    real(ReKi)              :: E          !< Exponential poart for the mollifider
    real(ReKi)              :: r3_inv     !< 
+   real(ReKi)              :: r2, r3, rc3!< |r|^2, |r|^3, RegParam^3 (reused to avoid recomputing ** intrinsics)
    real(ReKi)              :: rDeltaP    !< norm , distance between point and particle
    real(ReKi)              :: ScalarPart !< the part containing the inverse of the distance, but not 4pi, Mollifier
-   rDeltaP=sqrt(DeltaP(1)**2+ DeltaP(2)**2+ DeltaP(3)**2)! norm
+   r2      = DeltaP(1)**2+ DeltaP(2)**2+ DeltaP(3)**2
+   rDeltaP = sqrt(r2)! norm
    if (rDeltaP<MINNORM) then !--- Exactly on the Singularity 
       Ui(1:3)  = 0.0_ReKi
       return
    else !--- Normal Procedure 
+      r3 = r2*rDeltaP ! |r|^3, reused below
       C(1) = Alpha(2) * DeltaP(3) - Alpha(3) * DeltaP(2)
       C(2) = Alpha(3) * DeltaP(1) - Alpha(1) * DeltaP(3)
       C(3) = Alpha(1) * DeltaP(2) - Alpha(2) * DeltaP(1)
       select case (RegFunction) !
       case (idRegNone) ! No mollification
-         r3_inv     = 1._ReKi/(rDeltaP**3)
+         r3_inv     = 1._ReKi/r3
          ScalarPart = r3_inv*fourpi_inv
       case (idRegExp) ! Exponential mollifier
-         r3_inv     = 1._ReKi/(rDeltaP**3)
-         E          = exp(-rDeltaP**3/RegParam**3)
-         ScalarPart = (1._ReKi-E)*r3_inv*fourpi_inv
+         rc3        = RegParam*RegParam*RegParam
+         r3_inv     = 1._ReKi/r3
+         if (r3 > PART_REG_CUT3*rc3) then ! r > 2*rc: mollifier -> 1 (skip exp), consistent with far-field multipole floor
+            ScalarPart = r3_inv*fourpi_inv
+         else
+            E          = exp(-r3/rc3)
+            ScalarPart = (1._ReKi-E)*r3_inv*fourpi_inv
+         endif
       case (idRegCompact) ! Compact support
-         r3_inv     = 1._ReKi/sqrt(RegParam**6+rDeltaP**6)
+         rc3        = RegParam*RegParam*RegParam
+         r3_inv     = 1._ReKi/sqrt(rc3*rc3+r3*r3)
          ScalarPart = r3_inv*fourpi_inv
       case default 
          print*,'[ERROR] Wrong regularization function for particles',RegFunction

modules/aerodyn/src/FVW_Subs.f90

@@ -1047,7 +1047,8 @@ subroutine FVW_InitRegularization(x, p, m, ErrStat, ErrMsg)
          write(*,'(A,1F8.4)')   'RegParam (Input      ) : ',p%WakeRegParam
       endif
 
-      if (p%RegDeterMethod==idRegDeterConstant) then
+      select case (p%RegDeterMethod)
+      case (idRegDeterConstant)
          ! Constant reg param throughout the wake
          if (p%WakeRegMethod==idRegAge) then ! NOTE: age method implies a division by rc
             p%WingRegParam=max(0.01_ReKi, p%WingRegParam)
@@ -1061,7 +1062,7 @@ subroutine FVW_InitRegularization(x, p, m, ErrStat, ErrMsg)
             x%W(iW)%Eps_NW(1:3,:,2) = p%WakeRegParam ! Second age is always WakeRegParam
          endif
 
-      else if (p%RegDeterMethod==idRegDeterAuto) then
+      case (idRegDeterAuto)
          ! TODO this is beta
          print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          print*,'!!! NOTE: using optimized wake regularization parameters is still a beta feature!'
@@ -1078,14 +1079,14 @@ subroutine FVW_InitRegularization(x, p, m, ErrStat, ErrMsg)
          write(*,'(A,1F8.4)')   'WakeRegParam      : ', p%WakeRegParam
          write(*,'(A,1F8.4)')   'WingRegParam      : ', p%WingRegParam
          write(*,'(A,1F9.4)')   'CoreSpreadEddyVisc: ', p%CoreSpreadEddyVisc
-      ! Set reg param on wing and first NW
-      ! NOTE: setting the same in all three directions for now, TODO!
-      x%W(iW)%Eps_NW(1:3,:,1) = p%WingRegParam ! First age is always WingRegParam (LL)
-      if (p%nNWMax>1) then
-         x%W(iW)%Eps_NW(1:3,:,2) = p%WakeRegParam ! Second age is always WakeRegParam
-      endif
+         ! Set reg param on wing and first NW
+         ! NOTE: setting the same in all three directions for now, TODO!
+         x%W(iW)%Eps_NW(1:3,:,1) = p%WingRegParam ! First age is always WingRegParam (LL)
+         if (p%nNWMax>1) then
+            x%W(iW)%Eps_NW(1:3,:,2) = p%WakeRegParam ! Second age is always WakeRegParam
+         endif
 
-      else if (p%RegDeterMethod==idRegDeterChord) then
+      case (idRegDeterChord)
          ! Using chord to scale the reg param
          do iSpan=1,p%W(iW)%nSpan
             x%W(iW)%Eps_NW(1:3, iSpan, 1) = p%WingRegParam * p%W(iW)%chord_CP(iSpan)
@@ -1094,7 +1095,7 @@ subroutine FVW_InitRegularization(x, p, m, ErrStat, ErrMsg)
             endif
          enddo
 
-      else if (p%RegDeterMethod==idRegDeterSpan) then
+      case (idRegDeterSpan)
          ! Using dr to scale the reg param
          do iSpan=1,p%W(iW)%nSpan
             ds = p%W(iW)%s_LL(iSpan+1)-p%W(iW)%s_LL(iSpan)
@@ -1103,10 +1104,11 @@ subroutine FVW_InitRegularization(x, p, m, ErrStat, ErrMsg)
                x%W(iW)%Eps_NW(1:3, iSpan, 2) = p%WakeRegParam * ds
             endif
          enddo
-      else ! Should never happen (caught earlier)
+
+      case default ! Should never happen (caught earlier)
          ErrStat = ErrID_Fatal
          ErrMsg ='Regularization determination method not implemented'
-      endif
+      end select
 
       if (iW==1) then
       call WrScr(' - OLAF regularization parameters (for wing 1):')
@@ -1223,9 +1225,9 @@ subroutine SegmentsToPartWrap(Sgmt, nSeg, PartPerSegment, RegFunction, Part, all
       Part%RegFunction = idRegExp ! TODO need to find a good equivalence and potentially adapt Epsilon in SegmentsToPart
    endif
    if (DEV_VERSION) then
-      call find_nan_2D(Part%P    , 'SegmentsToPartWrap Part%P')
-      call find_nan_2D(Part%Alpha, 'SegmentsToPartWrap Part%Alpha')
-      if (any(Part%RegParam(:)<-9999.99_ReKi)) then
+      call find_nan_2D(Part%P(:,1:nPart)    , 'SegmentsToPartWrap Part%P')
+      call find_nan_2D(Part%Alpha(:,1:nPart), 'SegmentsToPartWrap Part%Alpha')
+      if (any(Part%RegParam(1:nPart)<-9999.99_ReKi)) then ! Only the active particles are filled; the preallocated tail keeps its sentinel
          print*,'Error in Segment to part conversion'
          STOP
       endif
@@ -1261,18 +1263,19 @@ subroutine InducedVelocitiesAll_Init(p, x, m, Sgmt, Part, Tree, Panl, ErrStat, E
    Sgmt%nAct  = nSeg
    Sgmt%nActP = nSegP
 
-   ! --- Convert to particles if needed
-   if ((p%VelocityMethod(iVel)==idVelocityTreePart) .or. (p%VelocityMethod(iVel)==idVelocityPart)) then
+   select case (p%VelocityMethod(iVel))
+   case (idVelocityPart)
+      ! --- Convert segments to particles
       call SegmentsToPartWrap(Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunction, Part, allocPart=allocPart)
-   endif
-
-   ! --- Grow tree if needed
-   if (p%VelocityMethod(iVel)==idVelocityTreePart) then
+   case (idVelocityTreePart)
+      ! --- Convert segments to particles
+      call SegmentsToPartWrap(Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunction, Part, allocPart=allocPart)
+      ! --- Grow particle tree
       call grow_tree_part(Tree, Part%nAct, Part%P, Part%Alpha, Part%RegFunction, Part%RegParam, 0)
-
-   elseif (p%VelocityMethod(iVel)==idVelocityTreeSeg) then
+   case (idVelocityTreeSeg)
+      ! --- Grow segment tree
       call grow_tree_segment(Tree, nSeg, Sgmt%Points, Sgmt%Connct(:,1:nSeg), Sgmt%Gamma(1:nSeg), p%RegFunction, Sgmt%Epsilon(1:nSeg), 0)
-   endif
+   end select
 
    ! --- Src
    Panl%p_Src => p%SrcPnl
@@ -1298,20 +1301,26 @@ subroutine InducedVelocitiesAll_Calc(CPs, nCPs, Uind, p, Sgmt, Part, Tree, Panl,
    ErrStat= ErrID_None
    ErrMsg =''
 
-   if (p%VelocityMethod(iVel)==idVelocityBasic) then
+   select case (p%VelocityMethod(iVel))
+   case (idVelocityBasic)
       call ui_seg( 1, nCPs, CPs, 1, Sgmt%nAct, Sgmt%Points, Sgmt%Connct, Sgmt%Gamma, Sgmt%RegFunction, Sgmt%Epsilon, Uind)
 
-   elseif (p%VelocityMethod(iVel)==idVelocityTreePart) then
+   case (idVelocityTreePart)
       ! Tree has already been grown with InducedVelocitiesAll_Init
       !call print_tree(Tree)
-      call ui_tree_part(Tree, nCPs, CPs, p%TreeBranchFactor(iVel), Tree%DistanceDirect, Uind, ErrStat, ErrMsg)
+      call ui_tree_part(Tree, nCPs, CPs, p%TreeBranchFactor(iVel), 0.0_ReKi, Uind, ErrStat, ErrMsg)
 
-   elseif (p%VelocityMethod(iVel)==idVelocityPart) then
+   case (idVelocityPart)
       call ui_part_nograd(nCPs, CPs, Part%nAct, Part%P, Part%Alpha, Part%RegFunction, Part%RegParam, Uind)
 
-   elseif (p%VelocityMethod(iVel)==idVelocityTreeSeg) then
-      call ui_tree_segment(Tree, CPs, nCPs, p%TreeBranchFactor(iVel), Tree%DistanceDirect, Uind, ErrStat, ErrMsg)
-   endif
+   case (idVelocityTreeSeg)
+      call ui_tree_segment(Tree, CPs, nCPs, p%TreeBranchFactor(iVel), 0.0_ReKi, Uind, ErrStat, ErrMsg)
+
+   case default
+      ErrStat = ErrID_Fatal
+      ErrMsg  = 'Velocity method not implemented'
+
+   end select
 
    ! --- Src Panels
    if (associated(Panl%p_Src)) then
@@ -1338,19 +1347,25 @@ subroutine InducedVelocitiesAll_End(p, Tree, Part, Panl, ErrStat, ErrMsg, deallo
    ErrStat= ErrID_None
    ErrMsg =''
 
-   if (p%VelocityMethod(iVel)==idVelocityBasic) then
+   select case (p%VelocityMethod(iVel))
+   case (idVelocityBasic)
       ! Nothing
 
-   elseif (p%VelocityMethod(iVel)==idVelocityTreePart) then
+   case (idVelocityTreePart)
       if (deallocPart) deallocate(Part%P, Part%Alpha, Part%RegParam)
       call cut_tree(Tree)
 
-   elseif (p%VelocityMethod(iVel)==idVelocityPart) then
+   case (idVelocityPart)
       if (deallocPart) deallocate(Part%P, Part%Alpha, Part%RegParam)
 
-   elseif (p%VelocityMethod(iVel)==idVelocityTreeSeg) then
+   case (idVelocityTreeSeg)
       call cut_tree(Tree) ! We do not deallocate segment
-   endif
+
+   case default
+      ErrStat = ErrID_Fatal
+      ErrMsg  = 'Velocity method not implemented'
+
+   end select
 
    ! Src Panels (we nullify only)
    nullify(Panl%p_Src)
@@ -1547,26 +1562,31 @@ subroutine LiftingLineInducedVelocities(p, x, InductionAtCP, iDepthStart, m, Err
 
       ! --- Compute velocity on LL
       ! TreeSeg is faster but introduce some noise, so we keep this open for the user to choose
-      if (p%VelocityMethod(iVel) == idVelocityBasic) then
+      select case (p%VelocityMethod(iVel))
+      case (idVelocityBasic)
          call ui_seg( 1, nCPs, CPs, 1, nSeg, m%Sgmt%Points, m%Sgmt%Connct, m%Sgmt%Gamma, m%Sgmt%RegFunction, m%Sgmt%Epsilon, Uind)
 
-      else if (p%VelocityMethod(iVel) == idVelocityPart) then
+      case (idVelocityPart)
          call SegmentsToPartWrap(m%Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunction, m%Part, allocPart=.false.)
          call ui_part_nograd(nCPs, CPs, m%Part%nAct, m%Part%P, m%Part%Alpha, m%Part%RegFunction, m%Part%RegParam, Uind)
          !deallocate(Part%P, Part%Alpha, Part%RegParam)
 
-      else if (p%VelocityMethod(iVel) == idVelocityTreeSeg) then
+      case (idVelocityTreeSeg)
          call grow_tree_segment(Tree, nSeg, m%Sgmt%Points, m%Sgmt%Connct(:,1:nSeg), m%Sgmt%Gamma(1:nSeg), m%Sgmt%RegFunction, m%Sgmt%Epsilon(1:nSeg), 0)
          call ui_tree_segment(Tree, CPs, nCPs, p%TreeBranchFactor(iVel), DistanceDirect, Uind, ErrStat, ErrMsg)
          call cut_tree(Tree)
 
-      else if (p%VelocityMethod(iVel) == idVelocityTreePart) then
+      case (idVelocityTreePart)
          call SegmentsToPartWrap(m%Sgmt, nSeg, p%PartPerSegment(iVel), p%RegFunction, m%Part, allocPart=.false.)
          call grow_tree_part(Tree, m%Part%nAct, m%Part%P, m%Part%Alpha, m%Part%RegFunction, m%Part%RegParam, 0)
          call ui_tree_part(Tree, nCPs, CPs, p%TreeBranchFactor(iVel), DistanceDirect, Uind, ErrStat, ErrMsg)
          !deallocate(Part%P, Part%Alpha, Part%RegParam)
          call cut_tree(Tree)
-      endif
+
+      case default
+         ErrStat = ErrID_Fatal
+         ErrMsg  = 'Velocity method not implemented'
+      end select
       ! --- Src Panel contribution
       if (p%SrcPnl%n>0) then
          call ui_quad_src_nn(CPs, m%SrcPnl%RHS, p%SrcPnl%xi, p%SrcPnl%eta, p%SrcPnl%Pcent, p%SrcPnl%R_g2p, Uind, nCPs, p%SrcPnl%n)

modules/aerodyn/src/FVW_VortexTools.f90

@@ -63,6 +63,7 @@ module FVW_VortexTools
       real(ReKi)                         :: radius !< Typical dimension of a cell (max of x,y,z extent)
       real(ReKi),dimension(3)            :: center !< Used to store first the geometric center, then the vorticity center
       real(ReKi),dimension(3,10)         :: Moments
+      real(ReKi)                         :: maxRegParam = 0.0_ReKi !< Max reg. param (eps) over all particles in this node
       integer,dimension(:),pointer       :: iPart=>null()  !< indexes of particles stored in this node
       integer,dimension(:),pointer       :: leaves=>null()  ! NOTE: leaves are introduced to save memory
       type(T_Node),dimension(:), pointer :: branches =>null()
@@ -71,11 +72,10 @@ module FVW_VortexTools
 
    !> The type tree contains some basic data, a chained-list of nodes, and a pointer to the Particle data that were used
    type T_Tree
-      type(T_VPart)  :: Part            !< Storage for all particles
-      type(T_VSgmt)   :: Seg             !< Storage for all segments
+      type(T_VPart) :: Part            !< Storage for all particles
+      type(T_VSgmt) :: Seg             !< Storage for all segments
       integer       :: iStep =-1       !< Time step at which the tree was built
       logical       :: bGrown =.false. !< Is the tree build
-      real(ReKi)    :: DistanceDirect
       type(T_Node)  :: Root            !< Contains the chained-list of nodes
    end type T_Tree
 
@@ -587,6 +587,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%maxRegParam = Part%RegParam(1)
          nullify(node%iPart)
          nullify(node%branches)
          allocate(node%leaves(1:1))
@@ -625,7 +626,6 @@ subroutine grow_tree_part(Tree, nPart, PartP, PartAlpha, PartRegFunction, PartRe
       endif
       Tree%iStep  = iStep
       Tree%bGrown = .true.
-      Tree%DistanceDirect = 2*sum(PartRegParam)/size(PartRegParam) ! 2*mean(eps), below that distance eps has a strong effect
    end subroutine grow_tree_part
 
    !> Recursive function to grow/setup a tree. 
@@ -675,13 +675,15 @@ subroutine grow_tree_part_substep(node, Part)
       GeomC = node%center ! NOTE: we rely on the fact that our parent has set this to the Geometric value 
       VortC = 0.0_ReKi
       wTot = 0.0_ReKi
+      node%maxRegParam = 0.0_ReKi
       ! --- Barycenter of vorticity of the node 
       do i = 1,node%nPart
          PartPos   = Part%P(:,node%iPart(i))
          PartAlpha = Part%Alpha(:,node%iPart(i))
          wLoc      = (PartAlpha(1)**2 + PartAlpha(2)**2 + PartAlpha(3)**2)**0.5_ReKi ! Vorticity norm
          VortC     = VortC + wLoc*PartPos                                   ! Sum coordinates weighted by vorticity
          wTot      = wTot + wLoc                                                     ! Total vorticity
+         node%maxRegParam = max(node%maxRegParam, Part%RegParam(node%iPart(i))) ! Regularization floor uses max eps in this cell
       end do
       ! There is no vorticity, we make it a empty node and we exit 
       if(EqualRealNos(abs(wTot),0.0_ReKi)) then
@@ -925,6 +927,7 @@ subroutine grow_tree_segment(Tree_Seg, nSeg, SegPoints, SegConnct, SegGamma, Seg
          node%radius=0
          node%center(1:3)=0.5_ReKi*(Seg%SP(1:3,1)+Seg%SP(1:3,2))
          node%Moments=0.0_ReKi
+         node%maxRegParam = Seg%RegParam(1)
          nullify(node%iPart)
          nullify(node%branches)
          allocate(node%leaves(1:1))
@@ -973,7 +976,6 @@ subroutine grow_tree_segment(Tree_Seg, nSeg, SegPoints, SegConnct, SegGamma, Seg
       endif
       Tree_Seg%iStep  = iStep
       Tree_Seg%bGrown = .true.
-      Tree_Seg%DistanceDirect = 2*sum(SegRegParam)/size(SegRegParam) ! 2*mean(eps), below that distance eps has a strong effect ! TODO REMOVE
    end subroutine grow_tree_segment
 
    !> Recursive function to grow/setup a tree.
@@ -1024,13 +1026,15 @@ subroutine grow_tree_segment_substep(node, Seg)
       GeomC = node%center ! NOTE: we rely on the fact that our parent has set this to the Geometric value
       VortC = 0.0_ReKi
       wTot = 0.0_ReKi
+      node%maxRegParam = 0.0_ReKi
       ! --- Barycenter of vorticity of the node
       do i = 1,node%nPart
          P1 = Seg%SP(1:3,Seg%SConnct(1,node%iPart(i)))
          P2 = Seg%SP(1:3,Seg%SConnct(2,node%iPart(i)))
          SegCenter = 0.5_ReKi*(P1+P2)
          VortC = VortC + abs(Seg%SGamma(node%iPart(i)))*SegCenter  ! Sum coordinates weighted by vorticity
          wTot   = wTot + abs(Seg%SGamma(node%iPart(i)))            ! Total vorticity
+         node%maxRegParam = max(node%maxRegParam, Seg%RegParam(node%iPart(i))) ! Regularization floor uses max eps in this cell
       end do
       ! There is no vorticity, we make it a empty node and we exit
       if(EqualRealNos(abs(wTot),0.0_ReKi)) then
@@ -1501,7 +1505,7 @@ recursive subroutine ui_tree_part_11(node, CP, Uind)
                enddo
             endif
          else
-            distDirect = max(BranchFactor*node%radius, DistanceDirect) ! Under this distance-> Direct eval., Above it -> quadrupole calculation
+            distDirect = max(BranchFactor*node%radius + 2.0_ReKi*node%maxRegParam, DistanceDirect) ! Direct eval. below this distance; floor = branch radius + 2*max(eps) in cell
             DeltaP  = - node%center + CP(1:3)                          ! Vector between the control point and the center of the branch
             r       = sqrt( DeltaP(1)**2 + DeltaP(2)**2 + DeltaP(3)**2)
             ! Test if the control point is too close from the branch node so that a direct evaluation is needed 
@@ -1704,7 +1708,7 @@ recursive subroutine ui_tree_segment_11(node, CP, Uind)
                enddo
             endif
          else
-            distDirect = max(BranchFactor*node%radius, DistanceDirect) ! Under this distance-> Direct eval., Above it -> quadrupole calculation
+            distDirect = max(BranchFactor*node%radius + 2.0_ReKi*node%maxRegParam, DistanceDirect) ! Direct eval. below this distance; floor = branch radius + 2*max(eps) in cell
             DeltaP  = - node%center + CP(1:3)                          ! Vector between the control point and the center of the branch
             r       = sqrt( DeltaP(1)**2 + DeltaP(2)**2 + DeltaP(3)**2)
             ! Test if the control point is too close from the branch node so that a direct evaluation is needed