|
| 1 | +# Calling CCBlade.jl from Python (optionally with derivatives) |
| 2 | + |
| 3 | +In this example we repeat the [quick start](tutorial.md), but with Python. We won't repeat all the usage details described in the quick start. We are using the [JuliaCall](https://juliapy.github.io/PythonCall.jl/stable/) packages to enable this. Notice that the main change from the Julia example is that we need to explicitly use the Julia broadcast functions (which is what the dot syntax does in Julia). |
| 4 | + |
| 5 | +```python |
| 6 | +import numpy as np |
| 7 | +from juliacall import Main as jl |
| 8 | +import matplotlib.pyplot as plt |
| 9 | + |
| 10 | +jl.seval("using CCBlade") |
| 11 | + |
| 12 | +Rtip = 10/2.0 * 0.0254 # inches to meters |
| 13 | +Rhub = 0.10*Rtip |
| 14 | +B = 2 # number of blades |
| 15 | + |
| 16 | +rotor = jl.Rotor(Rhub, Rtip, B) |
| 17 | + |
| 18 | +propgeom = np.array([ |
| 19 | + [0.15, 0.130, 32.76], |
| 20 | + [0.20, 0.149, 37.19], |
| 21 | + [0.25, 0.173, 33.54], |
| 22 | + [0.30, 0.189, 29.25], |
| 23 | + [0.35, 0.197, 25.64], |
| 24 | + [0.40, 0.201, 22.54], |
| 25 | + [0.45, 0.200, 20.27], |
| 26 | + [0.50, 0.194, 18.46], |
| 27 | + [0.55, 0.186, 17.05], |
| 28 | + [0.60, 0.174, 15.97], |
| 29 | + [0.65, 0.160, 14.87], |
| 30 | + [0.70, 0.145, 14.09], |
| 31 | + [0.75, 0.128, 13.39], |
| 32 | + [0.80, 0.112, 12.84], |
| 33 | + [0.85, 0.096, 12.25], |
| 34 | + [0.90, 0.081, 11.37], |
| 35 | + [0.95, 0.061, 10.19], |
| 36 | + [1.00, 0.041, 8.99]] |
| 37 | +) |
| 38 | + |
| 39 | +r = propgeom[:, 0] * Rtip |
| 40 | +chord = propgeom[:, 1] * Rtip |
| 41 | +theta = propgeom[:, 2] * np.pi/180 |
| 42 | + |
| 43 | +af = jl.AlphaAF("data/naca4412.dat") |
| 44 | + |
| 45 | +sections = jl.broadcast(jl.Section, r, chord, theta, jl.Ref(af)) |
| 46 | + |
| 47 | +Vinf = 5.0 |
| 48 | +Omega = 5400*np.pi/30 # convert to rad/s |
| 49 | +rho = 1.225 |
| 50 | + |
| 51 | +op = jl.broadcast(jl.simple_op, Vinf, Omega, r, rho) |
| 52 | + |
| 53 | +out = jl.broadcast(jl.solve, jl.Ref(rotor), sections, op) |
| 54 | + |
| 55 | +plt.figure() |
| 56 | +plt.plot(r/Rtip, out.Np) |
| 57 | +plt.plot(r/Rtip, out.Tp) |
| 58 | +plt.xlabel("r/Rtip") |
| 59 | +plt.ylabel("distributed loads (N/m)") |
| 60 | +plt.legend(["flapwise", "lead-lag"]) |
| 61 | +plt.show() |
| 62 | +``` |
| 63 | + |
| 64 | +We now continue the example demonstrating how to get derivatives for use in Python (but where the derivative computationa occurs in Julia via algorithmic differentiation). For this functionality we need to load the [PythonCall](https://juliapy.github.io/PythonCall.jl/stable/) package (which enables us to call back into Python from Julia), and we need the [ImplicitAD](https://github.com/byuflowlab/ImplicitAD.jl) packages which provides a convenience function for the derivative computation. Alternatively, for more advanced users you can just use the Julia differentiation packages directly (ForwardDiff, ReverseDiff, etc.) as demonstrated in the [how to guide](howto.md). Note that for Julia AD to work, all the function calls will need to be Julia function calls. Even though all the setup is happening here in Python, we are only setting up inputs. All functions are calls to Julia (jl.somefunction()) |
| 65 | + |
| 66 | +Not counting the one-time loading cost of starting up the Julia runtime, we find the computational time of computing derivatives to be essentially identical to computing in pure Julia (in other words there is minimal overhead between Python/Julia for these examples). The below example demonstrates a forward mode Jacobian and a forward mode Jacobian-vector product. Other options (reverse Jacobian, reverse vector-Jacobian product) are discussed in [ImplicitAD](https://github.com/byuflowlab/ImplicitAD.jl) docs. |
| 67 | + |
| 68 | +```python |
| 69 | +jl.seval("using PythonCall") |
| 70 | +jl.seval("using ImplicitAD: derivativesetup") |
| 71 | + |
| 72 | +nc = len(chord) |
| 73 | + |
| 74 | +# ImplicitAD expects a funciton we want to differentiate in the form f = func(x, p) |
| 75 | +# where f is output vector, x is input vector, and p are parameters we do not differentiate w.r.t. |
| 76 | +def ccbladewrap(x, p): |
| 77 | + chord = x[:nc] |
| 78 | + theta = x[nc:] |
| 79 | + sections = jl.broadcast(jl.Section, r, chord, theta, jl.Ref(af)) |
| 80 | + out = jl.broadcast(jl.solve, jl.Ref(rotor), sections, op) |
| 81 | + T, Q = jl.thrusttorque(rotor, sections, out) |
| 82 | + return [T, Q] |
| 83 | + |
| 84 | +x = np.concatenate([chord, theta]) |
| 85 | +p = () |
| 86 | +jacobian = jl.derivativesetup(ccbladewrap, x, p, "fjacobian") # a forward-mode Jacobian is one option |
| 87 | + |
| 88 | +# preallocate Jacobian then evaluate |
| 89 | +J = np.zeros((2, len(x))) |
| 90 | +jacobian(J, x) |
| 91 | +print(J) |
| 92 | +# can now change x, and evaluate jacobian(J, x) repeatedly at other points |
| 93 | + |
| 94 | + |
| 95 | +# demonstrate a Jacobian-vector product |
| 96 | +jvp = jl.derivativesetup(ccbladewrap, x, p, "jvp") |
| 97 | +xdot = np.ones(len(x)) |
| 98 | +fdot = np.zeros(2) |
| 99 | +jvp(fdot, x, xdot) |
| 100 | +print(fdot) |
| 101 | +# can continue to call jvp for different x, xdot pairs |
| 102 | +``` |
0 commit comments