Add a compact-support core regularization function for OLAF particle methods

docs/source/user/aerodyn-olaf/ExampleFiles/ExampleFile--OLAF.dat

@@ -22,7 +22,8 @@ default FWShedVorticity    - Include shed vorticity in the far wake {default: fa
 ------------------- WAKE REGULARIZATIONS AND DIFFUSION -----------------------------------------
 default DiffusionMethod    - Diffusion method to account for viscous effects {0: None, 1: Core Spreading, "default": 0}         
 2       RegDeterMethod     - Method to determine the regularization parameters {0: Constant, 1: Optimized, 2: Chord-scaled, 3: dr-scaled, default: 0 }          
-default RegFunction        - Viscous diffusion function {0: None, 1: Rankine, 2: LambOseen, 3: Vatistas, 4: Denominator, "default": 3} (switch)
+default RegFunction        - Segment regularization function {0: None, 1: Rankine, 2: LambOseen, 3: Vatistas, 4: Denominator, "default": 3} (switch)
+default RegFunctionPart    - Particle regularization function {0: None, 1: Exponential, 2: Compact, "default": 1} [only if VelocityMethod=2,3] (switch)
 default WakeRegMethod      - Wake regularization method {1: Constant, 2: Stretching, 3: Age, default: 1} (switch)
 0.25    WakeRegFactor      - Wake regularization factor (m or -)
 0.25    WingRegFactor      - Wing regularization factor (m or -)

docs/source/user/aerodyn-olaf/InputFiles.rst

@@ -178,13 +178,29 @@ See :numref:`Guidelines-OLAF` for recommendations on setting up this parameter.
 
 
 **RegFunction** [switch] specifies the regularization function used to remove
-the singularity of the vortex elements, as specified in
+the singularity of the vortex elements/segments, as specified in
 :numref:`sec:vortconv`. There are five options: 1) no correction *[0]*,
 2) the Rankine method *[1]*, 3) the Lamb-Oseen method *[2]*, 4) the Vatistas
 method *[3]*, and 5) the denominator offset method *[4]*. 
 The functions are given in :numref:`sec:RegularizationFunction`. 
 The default option is *[3]*.
 
+**RegFunctionPart** [switch] specifies the regularization function used for the
+vortex particles, which are used when a particle-based velocity method is
+selected (*VelocityMethod* = *[2,3]*). There are three options: 1) no
+correction *[0]*, 2) the exponential method *[1]*, and 3) the compact-support
+method *[2]*. The functions are given in
+:numref:`sec:RegularizationFunctionPart`. 
+The compact-support option *[2]* has a finite support radius beyond which the
+kernel reverts to the exact singular kernel; with the tree-accelerated particle
+method (*VelocityMethod* = *[2]*) this allows more particle interactions to be
+handled by the far-field multipole approximation, so it can be significantly
+faster than the exponential option depending on the case (wake size, particle
+count, and regularization parameter).
+The compact-support kernel is also purely polynomial and avoids evaluating a
+transcendental (exponential) function, which can further reduce cost.
+The default option is *[1]*.
+
 **WakeRegMethod** [switch] specifies the method of determining viscous core
 radius (i.e., the regularization parameter). There are three options: 1)
 constant *[1]*, 2) stretching *[2]*, and 3) age *[3]*. The methods are

docs/source/user/aerodyn-olaf/OLAFTheory.rst

@@ -489,8 +489,16 @@ refinement of this option will be considered in the future.
 
 .. _sec:RegularizationFunction:
 
-Implemented regularization functions
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Segment regularization functions
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The regularization functions described in this section apply to the vortex
+*segments* used to represent both the bound (blade) vorticity and the wake
+vorticity, i.e., they regularize the segment Biot-Savart kernel of
+Eq. :eq:`eq:BiotSavartSegment`. They are selected with the input
+**RegFunction**. The corresponding regularization functions used by the
+vortex-particle representation of the wake are described separately in
+:numref:`sec:RegularizationFunctionPart`.
 
 Several regularization functions have been
 developed (:cite:`olaf-Rankine58_1,olaf-Scully75_1,olaf-Vatistas91_1`).  At present, five
@@ -562,7 +570,108 @@ Here, the singularity is removed by introducing an additive factor in the
 denominator of Eq. :eq:`eq:BiotSavartSegment`, proportional to the filament
 length :math:`r_0`. In this case, :math:`F_\nu=1`. This method is found in the
 work of van Garrel (:cite:`olaf-Garrel03_1`).
+.. _sec:RegularizationFunctionPart:
+
+Particle regularization functions
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The regularization functions of :numref:`sec:RegularizationFunction` apply to
+the vortex segments. When a particle-based velocity method is selected
+(**VelocityMethod=[3]**, the direct vortex-particle method, or
+**VelocityMethod=[2]**, the tree-accelerated particle method), the wake segments
+are converted into vortex particles and a separate set of regularization
+functions applies. The particle regularization function is selected with the
+input **RegFunctionPart**.
 
+Segment-to-particle conversion
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Each wake segment of circulation :math:`\Gamma` and length vector
+:math:`\vec{l}=\vec{x}_2-\vec{x}_1` is divided into :math:`n_p`
+(**PartPerSegment**) equally spaced particles. The intensity of a particle
+:math:`\vec{\alpha}` (the vorticity integrated over the volume represented by
+the particle, :math:`\vec{\alpha}=\vec{\omega}\,dV`) is obtained from the parent
+segment as
+
+.. math::
+   \vec{\alpha} = \frac{\Gamma\,\vec{l}}{n_p}
+   :label: eq:SegToPart
+
+so that the :math:`n_p` particles of a segment sum to the total segment
+vorticity :math:`\Gamma\,\vec{l}`. The particles are placed at the centers of
+the :math:`n_p` sub-segments, and each particle inherits the regularization
+parameter (core radius :math:`r_c`) of its parent segment.
+
+Particle velocity kernel
+^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The velocity induced at a point :math:`\vec{x}` by a vortex particle of
+intensity :math:`\vec{\alpha}` located at :math:`\vec{x}_p` is the regularized
+point-vortex kernel
+
+.. math::
+   \vec{v}(\vec{x}) = \frac{1}{4\pi}\, g(\bar{r})\,
+   \frac{\vec{\alpha}\times\vec{r}}{r^3}
+   ,\qquad \vec{r}=\vec{x}-\vec{x}_p,\quad r=|\vec{r}|,\quad \bar{r}=\frac{r}{r_c}
+   :label: eq:BiotSavartParticle
+
+where :math:`g` is the particle regularization factor (the particle-method
+analog of :math:`F_\nu`) and :math:`r_c` is the particle core radius. Away from
+the particle center the induced velocity decays as :math:`1/r^2`. If no
+correction is used (**RegFunctionPart=[0]**), :math:`g=1` and the singular
+point-vortex kernel is recovered.
+
+Exponential
+^^^^^^^^^^^
+
+If the exponential method is used (**RegFunctionPart=[1]**, the default), the
+regularization factor is
+
+.. math::
+   g(\bar{r}) = 1 - \exp\!\left(-\bar{r}^{\,3}\right)
+   :label: eq:PartExp
+
+For :math:`\bar{r} > 2` (beyond two core radii) the exponential term is
+negligible and :math:`g` is set to :math:`1`.
+
+Compact support
+^^^^^^^^^^^^^^^^
+
+If the compact-support method is used (**RegFunctionPart=[2]**), the kernel is
+*exactly* singular (:math:`g=1`) beyond a finite support radius
+:math:`r_c'=1.6\,r_c`, and inside the support it is given by a polynomial
+mollifier
+
+.. math::
+   g(\bar{r}') = \begin{cases}
+   \dfrac{35\,\bar{r}'^{3} - 42\,\bar{r}'^{5} + 15\,\bar{r}'^{7}}{8}
+   & 0 \le \bar{r}' < 1 \\[2mm]
+   1 & \bar{r}' \ge 1
+   \end{cases}
+   ,\qquad \bar{r}' = \frac{r}{r_c'} = \frac{r}{1.6\,r_c}
+   :label: eq:PartCompact
+
+The polynomial is smooth (:math:`C^\infty`) inside the support; the overall
+regularization factor is globally :math:`C^2`, its regularity being limited by
+the match to the singular kernel at the support boundary, where :math:`g=1`,
+:math:`g'=0`, and :math:`g''=0` but :math:`g'''\neq 0` at :math:`\bar{r}'=1`. At
+the center the mollifier vanishes as :math:`\bar{r}'^{3}`
+(:math:`g=g'=g''=0`), so the induced velocity is regular there. The
+support-radius factor :math:`1.6` makes the compact kernel roughly equivalent
+to the exponential kernel with the same :math:`r_c`. It matches the peak
+tangential (swirl) velocity *magnitude* of the two kernels, while being a
+reasonable compromise across three possible equivalence criteria:
+
+.. math::
+   \frac{r_c'}{r_c} \approx \begin{cases}
+   1.54 & \text{radius of peak tangential velocity} \\
+   1.60 & \text{peak tangential velocity magnitude} \\
+   1.65 & \text{second moment of the vorticity distribution}
+   \end{cases}
+
+Because the kernel is exactly equal to the singular kernel beyond :math:`r_c'`,
+the tree-accelerated particle method can apply an exact far-field cutoff at
+that radius.
 .. _sec:corerad:
 
 Time Evolution of the Regularization Parameter–Core Spreading Method

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.f90

@@ -456,6 +456,7 @@ SUBROUTINE FVW_SetParametersFromInputFile( InputFileData, p, ErrStat, ErrMsg )
    p%FWShedVorticity      = InputFileData%FWShedVorticity
    p%DiffusionMethod      = InputFileData%DiffusionMethod
    p%RegFunction          = InputFileData%RegFunction
+   p%RegFunctionPart      = InputFileData%RegFunctionPart
    p%RegDeterMethod       = InputFileData%RegDeterMethod
    p%WakeRegMethod        = InputFileData%WakeRegMethod
    p%WakeRegParam         = InputFileData%WakeRegParam

modules/aerodyn/src/FVW_BiotSavart.f90

@@ -10,6 +10,9 @@ 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 :: PART_REG_C2   =  1.6_ReKi             !< Compact-C2 support = PART_REG_C2*RegParam (exp-equivalent core); also the exact multipole floor for that kernel
    real(ReKi),parameter :: MINDENOM=0.0_ReKi
 !    real(ReKi),parameter :: MINDENOM=1e-15_ReKi
    real(ReKi),parameter :: MINNORM=1e-4
@@ -29,6 +32,21 @@ module FVW_BiotSavart
 
 contains
 
+!> Multipole floor factor for the particle kernels: number of core radii beyond which the
+!! regularized kernel equals the singular 1/r^3, so the far-field multipole expansion is valid.
+pure function PartRegFloorFactor(RegFunction) result(f)
+   integer(IntKi), intent(in) :: RegFunction
+   real(ReKi)                 :: f
+   select case (RegFunction)
+   case (idRegNone)    ! Unregularized kernel: no extra core-based floor needed
+      f = 0.0_ReKi
+   case (idRegCompact) ! Truly compact: exactly singular beyond rc = PART_REG_C2*RegParam
+      f = PART_REG_C2
+   case default        ! Exponential: conservative 2*rc floor
+      f = PART_REG_NRAD
+   end select
+end function PartRegFloorFactor
+
 
 !> Induced velocity from one segment at one control points
 subroutine ui_seg_11(DeltaPa, DeltaPb, SegGamma, RegFunction, RegParam1, Uind)
@@ -374,27 +392,43 @@ 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, core^3 (reused to avoid recomputing ** intrinsics)
+   real(ReKi)              :: rc2, t     !< compact-C2 core^2 and (r/rc)^2
    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
-      case (idRegCompact) ! Compact support
-         r3_inv     = 1._ReKi/sqrt(RegParam**6+rDeltaP**6)
-         ScalarPart = 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) ! Truly compact C2 core: exactly singular for r>=rc, rc=PART_REG_C2*RegParam
+         rc2 = (PART_REG_C2*RegParam)**2
+         if (r2 >= rc2) then ! outside core: kernel is exactly the singular 1/r^3
+            r3_inv     = 1._ReKi/r3
+            ScalarPart = r3_inv*fourpi_inv
+         else                ! inside core: C2 polynomial mollifier
+            t          = r2/rc2
+            rc3        = rc2*PART_REG_C2*RegParam
+            ScalarPart = (35._ReKi + t*(-42._ReKi + 15._ReKi*t))*fourpi_inv/(8._ReKi*rc3)
+         endif
       case default 
          print*,'[ERROR] Wrong regularization function for particles',RegFunction
          STOP

modules/aerodyn/src/FVW_IO.f90

@@ -60,6 +60,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg )
    CALL ReadVarWDefault(UnIn,FileName,Inp%DiffusionMethod    ,'DiffusionMethod'   ,'',idDiffusionNone , ErrStat2,ErrMsg2); if(Failed())return
    CALL ReadVarWDefault(UnIn,FileName,Inp%RegDeterMethod     ,'RegDeterMethod'    ,'',idRegDeterConstant, ErrStat2,ErrMsg2); if(Failed())return
    CALL ReadVarWDefault(UnIn,FileName,Inp%RegFunction        ,'RegFunction'       ,'',idRegVatistas   , ErrStat2,ErrMsg2); if(Failed())return
+   CALL ReadVarWDefault(UnIn,FileName,Inp%RegFunctionPart    ,'RegFunctionPart'   ,'',idRegExp        , ErrStat2,ErrMsg2); if(Failed())return
    CALL ReadVarWDefault(UnIn,FileName,Inp%WakeRegMethod      ,'WakeRegMethod'     ,'',idRegAge        , ErrStat2,ErrMsg2); if(Failed())return
    CALL ReadVar        (UnIn,FileName,Inp%WakeRegParam       ,'WakeRegParam'      ,''                 , ErrStat2,ErrMsg2); if(Failed())return
    CALL ReadVar        (UnIn,FileName,Inp%WingRegParam       ,'WingRegParam'      ,''                 , ErrStat2,ErrMsg2); if(Failed())return
@@ -190,6 +191,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg )
    if (Check(.not.(ANY(idDiffusionVALID==Inp%DiffusionMethod)) , 'Diffusion method (DiffusionMethod) not implemented: '//trim(Num2LStr(Inp%DiffusionMethod)))) return
    if (Check(.not.(ANY(idRegDeterVALID ==Inp%RegDeterMethod))  , 'Regularization determination method (RegDeterMethod) not yet implemented: '//trim(Num2LStr(Inp%RegDeterMethod)))) return
    if (Check(.not.(ANY(idRegVALID      ==Inp%RegFunction  )), 'Regularization function (RegFunction) not implemented: '//trim(Num2LStr(Inp%RegFunction)))) return
+   if (Check(.not.(ANY(idRegPartVALID  ==Inp%RegFunctionPart)), 'Particle regularization function (RegFunctionPart) not implemented: '//trim(Num2LStr(Inp%RegFunctionPart)))) return
    if (Check(.not.(ANY(idRegMethodVALID==Inp%WakeRegMethod)), 'Wake regularization method (WakeRegMethod) not implemented: '//trim(Num2LStr(Inp%WakeRegMethod)))) return
    if (Check(.not.(ANY(idShearVALID    ==Inp%ShearModel   )), 'Shear model (ShearModel) not valid: '//trim(Num2LStr(Inp%ShearModel)))) return
    if (Check(.not.(ANY(idVelocityVALID ==Inp%VelocityMethod(1))), 'Velocity method (VelocityMethod(1)) not valid: '//trim(Num2LStr(Inp%VelocityMethod(1))))) return

modules/aerodyn/src/FVW_Registry.txt

@@ -116,6 +116,7 @@ typedef     ^                   ^               IntKi
 typedef     ^                   ^               ReKi                             CoreSpreadEddyVisc -  - -   "Eddy viscosity used in the core spreading method" 
 typedef     ^                   ^               IntKi                            RegDeterMethod   -  -  -  "Regularization determinatino method (manual, automatic)"          -
 typedef     ^                   ^               IntKi                            RegFunction   -  -  -  "Type of regularizaion function (LambOseen, Vatistas, see FVW_BiotSavart)"          -
+typedef     ^                   ^               IntKi                            RegFunctionPart -  -  -  "Type of particle regularization function (None, Exp, Compact; see FVW_BiotSavart)"          -
 typedef     ^                   ^               IntKi                            WakeRegMethod -  - -   "Method for regularization (constant, stretching, age, etc.)"  - 
 typedef     ^                   ^               ReKi                             WakeRegParam -  - -   "Initial value of the regularization parameter" 
 typedef     ^                   ^               ReKi                             WingRegParam -  - -   "Regularization parameter of the wing" 
@@ -334,6 +335,7 @@ typedef     ^                   ^               IntKi
 typedef     ^                   ^               ReKi                             CoreSpreadEddyVisc -  - -   "Eddy viscosity used in the core spreading method" 
 typedef     ^                   ^               IntKi                            RegDeterMethod   -  -  -  "Regularization determinatino method (manual, automatic)"          -
 typedef     ^                   ^               IntKi                            RegFunction   -  -  -  "Type of regularizaion function (LambOseen, Vatistas, see FVW_BiotSavart)"          -
+typedef     ^                   ^               IntKi                            RegFunctionPart -  -  -  "Type of particle regularization function (None, Exp, Compact; see FVW_BiotSavart)"          -
 typedef     ^                   ^               IntKi                            WakeRegMethod -  - -   "Method for regularization (constant, stretching, age, etc.)"  - 
 typedef     ^                   ^               ReKi                             WakeRegParam  -  - -   "Factor used in the regularization " 
 typedef     ^                   ^               ReKi                             WingRegParam  -  - -   "Factor used in the regularization "