Add improved actuator zone model

[anilyil]

Purpose

This adds the improved actuator zone model that includes more parameters to control the distribution of the actuator zone terms, as well as adding a swirl term. I believe some of the changes originated from the work of @shamsheersc19.

It is WIP right now but opening the PR to start the discussions on this.

Expected time until merged

Type of change

• [ ] Bugfix (non-breaking change which fixes an issue)
• [ ] New feature (non-breaking change which adds functionality)
• [ ] Breaking change (non-backwards-compatible fix or feature)
• [ ] Code style update (formatting, renaming)
• [ ] Refactoring (no functional changes, no API changes)
• [ ] Documentation update
• [ ] Maintenance update
• [ ] Other (please describe)

Testing

Checklist

• [ ] I have run flake8 and black to make sure the Python code adheres to PEP-8 and is consistently formatted
• [ ] I have formatted the Fortran code with fprettify or C/C++ code with clang-format as applicable
• [ ] I have run unit and regression tests which pass locally with my changes
• [ ] I have added new tests that prove my fix is effective or that my feature works
• [ ] I have added necessary documentation

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 50.00000% with 2 lines in your changes missing coverage. Please review.

Project coverage is 41.20%. Comparing base (8f6696b) to head (77745ae).

Files with missing lines	Patch %	Lines
adflow/pyADflow.py	50.00%	2 Missing :warning:

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #375   +/-   ##
=======================================
  Coverage   41.19%   41.20%           
=======================================
  Files          13       13           
  Lines        4127     4131    +4     
=======================================
+ Hits         1700     1702    +2     
- Misses       2427     2429    +2     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

:rocket: New features to boost your workflow:

• :snowflake: Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[eytanadler]

Currently, the actuator zone's thrust cannot be set through the AeroProblem as expected if the actuator zone uses the simpleProp type. The simpleProp uses the thrust attribute from actuatorRegionData rather than the force one in the residual routine. However, only the force attribute is set in the update BC call. To fix this, the thrust attribute must also be updated with something like actuatorRegions(iRegion)%thrust = bcDataIn(iVar). This could be done in an if/else based on the actuator region type. I'd make this change, but I don't have Tapenade set up to differentiate it. @eirikurj do you have an environment set up where you could easily make this change?

[eirikurj]

@eytanadler I have everything setup and can add this, but I dont have any case available to test this. Do you have anything that can be considered as a unit/regression test?

[eytanadler]

I don't have anything unit testable but the basic behavior necessary is that this wasn't working before:

solver = ADFLOW(options=options)
solver.addActuatorRegion(
    "az_surf.x",
    actType="simpleProp",
    axis1=np.array([0.0, 0.0, 0.0]),
    axis2=np.array([1.0, 0.0, 0.0]),
   familyName="az",
    thrust=0.0,
    propRadius=1.0,
    spinnerRadius=0.1,
)

# Doing this should make thrust settable through the AeroProblem
ap.setBCVar("Thrust", 0.0, "az")
ap.addDV("Thrust", family="az", name="az_thrust")
ap.setDesignVars({"az_thrust": 1e4})

solver(ap)

Specifically, adding the design variable to the AeroProblem and setting it should change the thrust (observable via flowpower cost function), but it doesn't without the changes I mentioned. I think something like this would make a fine unit test, right?

[shamsheersc19]

I'm not sure if you deleted the comment about the differences in the formulas in the code vs the paper, Eytan. I cannot find it on GitHub, but I see it here in my email. In fact = one / propRadius, I think I had the extra radius in the denominator to make the very small numbers larger for the summing and dividing later. In the end, the extra divide by propRadius (which is the outer radius i.e., a constant value) gets cancelled out in the step later where region%thrustVec = region%thrustVec / region%totalThrustSum. The (two * pi * region%cellRadii(region%nCellIDs)) is to convert the force per unit radius (which is what the formula gives) to per unit volume, as described in the paper. I haven't included the disk thickness in the denominator in the code because it is a constant that cancels out in the later step I mentioned earlier.

I need to think more about what consequences that dot product error that you fixed will have.

On Fri, Mar 7, 2025 at 1:49 PM Eytan Adler @.***> wrote:

I just noticed that I don't think the mDistribParam and mDistribParam match what is in Sham's paper https://www.researchgate.net/profile/Shamsheer-Chauhan/publication/350073795_RANS-Based_Aerodynamic_Shape_Optimization_of_a_Wing_Considering_Propeller-Wing_Interaction/links/61630b65ae47db4e57bc9760/RANS-Based-Aerodynamic-Shape-Optimization-of-a-Wing-Considering-Propeller-Wing-Interaction.pdf. His paper has essentially $f_x = r^m (1 - r)^n$ (I'm simplifying here), but the code https://github.com/anilyil/adflow/blob/4745f0cea2183c56496143c435abe4f0ee6bbace/src/solver/actuatorRegion.F90#L279 looks like it implements something more like $f_x = r^m (1 - r)^n / r$ with the extra region%cellRadii(region%nCellIDs) in the denominator. This reduces $m$ by 1 relative to the equation in the paper. This isn't necessarily a problem, but it's not intuitive. I would suggest either making the small tweak to the implementation to match the paper or changing the default mDistribParam to be 2.0 instead of 1.0 and putting a comment in the docstring because the default distribution right now not very reasonable. Am I missing something here?

— Reply to this email directly, view it on GitHub https://github.com/mdolab/adflow/pull/375#issuecomment-2707177036, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADLM7IPHVQDJEDRRF23DAC32THS4RAVCNFSM6AAAAABVK6YS3OVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDOMBXGE3TOMBTGY . You are receiving this because you were mentioned.Message ID: @.***> [image: eytanadler]eytanadler left a comment (mdolab/adflow#375) https://github.com/mdolab/adflow/pull/375#issuecomment-2707177036

I just noticed that I don't think the mDistribParam and mDistribParam match what is in Sham's paper https://www.researchgate.net/profile/Shamsheer-Chauhan/publication/350073795_RANS-Based_Aerodynamic_Shape_Optimization_of_a_Wing_Considering_Propeller-Wing_Interaction/links/61630b65ae47db4e57bc9760/RANS-Based-Aerodynamic-Shape-Optimization-of-a-Wing-Considering-Propeller-Wing-Interaction.pdf. His paper has essentially $f_x = r^m (1 - r)^n$ (I'm simplifying here), but the code https://github.com/anilyil/adflow/blob/4745f0cea2183c56496143c435abe4f0ee6bbace/src/solver/actuatorRegion.F90#L279 looks like it implements something more like $f_x = r^m (1 - r)^n / r$ with the extra region%cellRadii(region%nCellIDs) in the denominator. This reduces $m$ by 1 relative to the equation in the paper. This isn't necessarily a problem, but it's not intuitive. I would suggest either making the small tweak to the implementation to match the paper or changing the default mDistribParam to be 2.0 instead of 1.0 and putting a comment in the docstring because the default distribution right now not very reasonable. Am I missing something here?

— Reply to this email directly, view it on GitHub https://github.com/mdolab/adflow/pull/375#issuecomment-2707177036, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADLM7IPHVQDJEDRRF23DAC32THS4RAVCNFSM6AAAAABVK6YS3OVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDOMBXGE3TOMBTGY . You are receiving this because you were mentioned.Message ID: @.***>

[shamsheersc19]

I had a look at the dot product error fixed in commit 0d21897. I got lucky there because in the cases I ran for the papers, v1(1) = 0 (or approximately 0) so it being multiplied by the incorrect number fortunately didn't affect my cases. Nice catch!

[eytanadler]

@shamsheersc19 I deleted it because I realized I overlooked the force per length to force per volume conversion and it made sense after. Thanks for the clarification!