From 8338af108c3876c9d7adaa10dc41798bf6676a0d Mon Sep 17 00:00:00 2001 From: swryan <881430+swryan@users.noreply.github.com> Date: Thu, 31 Oct 2024 20:47:37 -0400 Subject: [PATCH] Changed 'latest' job in test and docs workflows to use Python 3.12 (#1122) * Changed 'latest' job in test workflow to use Python 3.12 * Changed 'latest' job in docs workflow to use Python 3.12 * om_version check, doc build with latest * update cannonball examples to use OpenMDAO InterpND * adjust tolerance on cannonball load_case test * update phase to use OM InterpND; add python 3.13 job to test workflow * revert phase change * fix for failing balanced field length tests. * more generalization of the balanced field example * PEP, require_pyoptsparse * require_pyoptsparse --------- Co-authored-by: swryan Co-authored-by: Rob Falck --- .github/workflows/dymos_docs_workflow.yml | 6 +- .github/workflows/dymos_tests_workflow.yml | 26 +- .../balanced_field_funccomp.ipynb | 6 +- .../brachistochrone/brachistochrone_fbd.png | Bin 12202 -> 12188 bytes .../cannonball_implicit_duration.ipynb | 24 +- .../balanced_field/balanced_field_length.py | 189 ++++++++ .../test/test_balanced_field_func_comp.py | 193 +------- .../test/test_balanced_field_length.py | 424 +----------------- dymos/examples/cannonball/cannonball_ode.py | 14 +- .../test_doc_cannonball_implicit_duration.py | 21 +- .../doc/test_doc_two_phase_cannonball.py | 18 +- .../test/test_load_case_missing_phase.py | 2 +- .../test_two_phase_cannonball_birkhoff.py | 18 +- pyproject.toml | 3 + 14 files changed, 304 insertions(+), 640 deletions(-) create mode 100644 dymos/examples/balanced_field/balanced_field_length.py diff --git a/.github/workflows/dymos_docs_workflow.yml b/.github/workflows/dymos_docs_workflow.yml index 07956a67a..7b89d6eb0 100644 --- a/.github/workflows/dymos_docs_workflow.yml +++ b/.github/workflows/dymos_docs_workflow.yml @@ -31,14 +31,15 @@ jobs: SCIPY: '1.13' PETSc: '3.19' PYOPTSPARSE: 'v2.11.0' - OPENMDAO: 'dev' + OPENMDAO: 'latest' OPTIONAL: '[docs]' JAX: '0.4.28' PUBLISH_DOCS: 1 # make sure the latest versions of things don't break the docs + # sticking with Python 3.12 for now, 3.13 requires NumPy 2.1 which does not work yet with PETSc/pyoptsparse - NAME: latest - PY: 3 + PY: '3.12' NUMPY: 1 SCIPY: 1 PETSc: 3 @@ -268,7 +269,6 @@ jobs: cd $HOME/work/dymos/dymos ghp-import -n -p -f docs/dymos_book/_build/html - - name: Publish docs to openmdao.org if: | github.event_name == 'push' && github.ref == 'refs/heads/master' && diff --git a/.github/workflows/dymos_tests_workflow.yml b/.github/workflows/dymos_tests_workflow.yml index 3a1b9cbde..e6b374099 100644 --- a/.github/workflows/dymos_tests_workflow.yml +++ b/.github/workflows/dymos_tests_workflow.yml @@ -60,6 +60,12 @@ on: required: false default: false + python313: + type: boolean + description: "Include 'python313' in test matrix" + required: false + default: false + oldest: type: boolean description: "Include 'oldest' in test matrix" @@ -134,8 +140,9 @@ jobs: EXCLUDE: ${{ github.event_name == 'workflow_dispatch' && ! inputs.no_mpi }} # try latest versions + # sticking with Python 3.12 here, 3.13 requires NumPy 2.1 which does not work yet with pyoptsparse - NAME: latest - PY: 3 + PY: '3.12' NUMPY: 1 SCIPY: 1 PETSc: 3.21.0 @@ -146,6 +153,19 @@ jobs: JAX: 'latest' EXCLUDE: ${{ github.event_name == 'workflow_dispatch' && ! inputs.latest }} + # Python 3.13 (requires NumPy 2.1 which does not work yet with pyoptsparse) + - NAME: python313 + PY: '3.13' + NUMPY: 2 + SCIPY: 1 + PETSc: 3 + # PYOPTSPARSE: 'latest' + SNOPT: 7.7 + OPENMDAO: 'dev' + OPTIONAL: '[test]' + JAX: 'latest' + EXCLUDE: ${{ github.event_name == 'workflow_dispatch' && ! inputs.python313 }} + # oldest supported versions - NAME: oldest PY: 3.9 @@ -237,9 +257,9 @@ jobs: if [[ "${{ matrix.OPENMPI }}" && "${{ matrix.MPI4PY }}" ]]; then conda install openmpi=${{ matrix.OPENMPI }} mpi4py=${{ matrix.MPI4PY }} petsc4py=${{ matrix.PETSc }} -q -y elif [[ "${{ matrix.MPI4PY }}" ]]; then - conda install mpi4py=${{ matrix.MPI4PY }} petsc4py=${{ matrix.PETSc }} -q -y + conda install mpich mpi4py=${{ matrix.MPI4PY }} petsc4py=${{ matrix.PETSc }} -q -y else - conda install mpi4py petsc4py=${{ matrix.PETSc }} -q -y + conda install mpich mpi4py petsc4py=${{ matrix.PETSc }} -q -y fi export OMPI_MCA_rmaps_base_oversubscribe=1 export OMPI_MCA_btl=^openib diff --git a/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb b/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb index 829c43ec1..c89dca539 100644 --- a/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb +++ b/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb @@ -168,6 +168,7 @@ "source": [ "import openmdao.api as om\n", "import openmdao.func_api as omf\n", + "from dymos.utils.misc import om_version", "\n", "def wrap_ode_func(num_nodes, flight_mode, grad_method='jax', jax_jit=True):\n", " \"\"\"\n", @@ -223,7 +224,10 @@ " meta.declare_coloring('*', method=grad_method)\n", " meta.declare_partials(of='*', wrt='*', method=grad_method)\n", " \n", - " return om.ExplicitFuncComp(meta, use_jax=grad_method == 'jax', use_jit=jax_jit)\n", + " if om_version()[0] > (3, 35, 0):\n", + " return om.ExplicitFuncComp(meta, derivs_method=grad_method, use_jit=jax_jit)\n", + " else:\n", + " return om.ExplicitFuncComp(meta, use_jax=grad_method == 'jax', use_jit=jax_jit)\n", " \n", " " ] diff --git a/docs/dymos_book/examples/brachistochrone/brachistochrone_fbd.png b/docs/dymos_book/examples/brachistochrone/brachistochrone_fbd.png index a431c9811f59969a4aee1a22da8f81a85aa367f1..711521e8dc380844ece271e9b7fdd310b55629d4 100644 GIT binary patch literal 12188 zcmeHtX*ksH|F0zp;VWBZs|+&MkY#K|ma!JHWSx+m!DL^veM4lJ#x|rZNw#F)O;Tiu z>`~yq!GUB}L^# zC4{dzqEJY0n3$N`|9U{w!^=VJx*p{yScC?tY35Bv#)u&Ok-t+ea3UiUqk-SMYvh+n zoD9r-W&5jRqgT>USL%6>nkHWg@n*~zQS*VIT-TCDG3JG32*YxN>B^@>`m!M(Xw9B z5rD@=8Lm<}cD)i&*!avm;oYnFJyT9{?Uu`su=CJ|Z<4(V$*6m914yQq*d4LvV~K+7 z5SVq)M!UZHFB9wF+p`^s4r8C3X!4RF?VnY|?pW6Q@X5$vnRwNYe|vj0?yzZa`$H<8 zbpJAGK72^A1a@8iH-Z=bD~=o7+}PbF>oA_PU&`GlVKCI<9UmW`{v-&TFO_%OQrw}3 zlR;_U2IjN!6)u3gSNy2bZ=;<(Ie`xjZEd*eE-*JgZ_pePtaP~emYCJs+4=bN_<$4H zLQcb%48hnX^!fHj(lCeIpLd{0q%8BCRmXZ^AbelE7v-I1+M=eo1}V6yIw|ns?51Tl zcOTb+w`TBe<4}Efse^sDX^;!P9Vg+NKMM^7l624ze}-TIw@?N}x>Wj~ypmfXY{H9w zb8jDn~1AS#T?u|5D5v{wKX^4dsUNc*yx`ugnv=+1CnyPSxJK!d1 zyCdM!T;pQ5y?hWF!sf4Bx#0#A3Ncn;Rqj!56P{I1MH}izEoAG1jR_jFf#{WEPM2Yv zQ$pQD5z|pO;0Ey^%xqq8c$mwT=x3J^>LKuhZl-|irqv7#$_f)1gle(lSsu6m8q)qn zMGR8fboJFU-Mlu;ng$P!^m85Xb4Uj_-2#q!NoxBKRfaSl=>*JljsHw~Ga40AY|Xdg zeEEU7SeU_KM=N^U(j21g5e_|(cm|TW1kZy^9bGC=UFboWjaW6d5na&`q9|SC%RB>d z8?qEJMvn&avtx(;=OB`ior>4b4oyZRgL@Sbs}{2BExF8abgU$u;`Q3vzABN1-A6SP zA=``M0a%EW32@KL_=|q~h@0lp6K|D}Hl#J$)1tCo2kcD0GH=_M!{QyX`7$)nMnQ6O zS!E*TZSOrjn&}0yE{HX+{H%vKEy$_pn-&_ZkJ!#MZRe=sga3fN3DY%Ra2)x^z2-XS z9+56PG>c^OjoFkyrgDEtPS11S)2RxwXnU2jN%rylb)oU!o0c%0Mu#VDbYd;jtt4Br z3a|w*_2<%9c1nONW6JZ*GHLX=hHZ8zFsOZJ4L?Qr`pbT=;?S_ho@OWGrvH!<>x zL!COvI&4SSyQ77-A474`l06sXKMoeLvt@J&m5<1vR;u+m zEGd)u@Js>chDsT$%ephDC+k|Z#cIRjKsib1<1OppUjJ!06S93>b$TWySr2dxSGEAi zRhek{)e9Y=`;_0av+b_J4G`X)iZ^O&krSTc?wnFCY^(ZWhWw=6wM}18ym93Uf)X3N z8@@C!lb?YmEg}ro-u@Mb>=c}7tJ#>%(aj)*J#9Mg_61y!goPi%RCjfy0@S|0uXF@Y zSNe!ZxzX|jkxka)o6m^Fu4ooXjm0FW2}71Y|3v5RX4iV?8J2YNyC3^_hdv2T_&dm+ z&hj8uCrCb^O)y+d>guLnu21Bqr{`!h_lSyegup#LBe`qWJnV}ev4Q{L9JKmS&7QZ$z#1D^+i?UMT^yRb&ek$9vu439C`$HC;8#FJseaUX?$^K3$;*XK+ zp<}A-Eek3!IP{b=2YW3~JH~tzY7{0z2;=b3A=QOuSrxJH^NgT=(RY+KP~c#sME!IG z9==nMNfCS)IqTgMo~NA*95oi4(5L+l?9+ccO>#!2nZ~^26%>Pl!uw0SBCS*~B6PjV zjH9FLV=l-p5mNpU88jC8CM0`X%p9U#n5(~_t*!=Tgon=30>z}Tf6b_bGk? zIQtxz`mLqnN3HuCgltiG`-_D>wb763ZFFVBR_lk$$djGFD++5Yc5eNOn|UQ22E6~7 znY6FNInErR-fseUdu-IKcW+b_-T(a5MLt9N?()$3d`4d`jZWI7{NH?VwG3XATY*j* zQmN6)KJ^ym;jS*Q!xFZvZ-@DaNUg8TwW=GskMLV5Cwc#rgU(VOwNk#mpd2%+6x`-$ zi8ey{GgzKIa2G`pX?ZnNh-2%!%$1jRDGLoAsv9{}bjr4%Elp*@^v_JYo;Num2J8>tcd+WWS z(R_~!bR!}d^+?+bWUC7)*Ufkxd{jm_Lc~P1J9J|8!OAoWF(tV5lQ)vLUD_J_kp%x4 zEkC%yrZ9?_@@EjwFKXWjfbTCfXrj;Q)w-}N9*b|Uav3LcNMD2jjt3yqm95Ta6d-8x zDXS2@-kZ-vj^sjVGwsRsQDHA$7m#QE!Plk@a?x8iuh5&-e%qA-9tL1 zeX<@{7v?uaP!O!6J@|fovx}6VokB=rC2(0T`kTtIyUQ30cd!)a8PNmrP!4V|^RW*j z)+1{>zSv+=TDBMuDO~q^kh!VrF1o!)t{M&E9N|Cms1mKm2j`23DhaB5LjTE?X03@% zH$75lx(lNXPTt-t=sQ=m;(8gz^4=)?iewvun<%3YoT^vFq)C1{tW4N7O#&gdvs_97 z<#zCCxZ|jVFe4qZgKK?%7L*tgs6)jUGAoFkdxS>(0bo0=>u+RP-%VJ^5Kl&y|dZRd~p<9j2{JC(z2&R#wAaD!8@jN?`5NOr5R?>c^)@S1~p5M%}l!YJ1F zvk4L8lx-U#)h)J(q#S&V$SLaES#2cxR85_-eTj!lCqmj=iY74sqOSHRab$&D?Zh(f zy_0cH0T~Ew5=QX>Ok4a(F-b!1 z&j1G_x(3;;Belhz&Rl_yU16(>jhA+_>AR+xHyutrJ!DR*kL@v0B{$|i4NNy3wQcCi zv^Jk&0Dwb9g|{b>YmIyT)gIQ;9)wI($~(aKCXlKCpapbV0SabS7ZcwsWCbQj_Q+C( zRrJ8(3bF#AGxt2<1v9?Q5!5^v7X~#FB6o=-zo=oOn&PhBcK%!sf}@2yaHVx2J5=vr zbBGIdkY#NL5f&u~y0`?4LkK1WAmZ%`|9P065=jo>_dR1yecg_B;=^5FSJXYz_0hoJL zpeRNwNw1*LwEPB7=<%)5(NUVbxN(Q<*3wXVI`-(NQv!OCotVGJi~!7C%56-FYQHo{>`t=jt%knvC@-!Frf zv*U$GW{j|;gUWan;%Lx^XtR$Tu}vfeiWg^2k*V{9>`LD&9Q|yEY9e>pe_~$_1Vr)M z1X`i5#9D_D!xIn$0!8cj1M7n(WqzLfi;*e0njrOtLQY&j@^-$W=jnwl*BhM8bL<~Y z*ZFB>&G%vkW35oCN`?%acOX%-%?)L{LyqP(j?|b?%p`#>3#(|`#iCHDf3oz16r-(Gph{F{w-ny~EN1O)^C--@TOOOO|4$7-Dn z8%$}PcQQZ`MzZhX^I7$j;Ltl@!~NyF@t2vTHT3gB>OSj0)nJ&%hEMnf#Q(iz@ws+^ z`>QJB3;2y<$Q(esjinJLr4yFJy@jI2-LX?|YOZ_lLqoMfo8l#KQ|H*oclL=MjB`9^ zWk2~?t-;ct_6W$Oya~TASBxka&dvY+-4esomWMTXZ@>qdMs7XWlGU{0Yc{OqP>TGj z-|f~u^rvNTLJ`bks`XoCPuCi&>(d*J3Jn_8+Zr=N4wb{5?%0BrL;(0K0R`HwYg zP1fG4wH^>IU3ADTI?Ed=hY^f~j&^qabEi~x;*NbkP^${~7b6f>I<~rPoDsBL`RssE zv~fFlTAfTi!BwTmc&*&u`{O8Nskyu_k4H1lV0$f8%*3LcVUga|zB5JD>F1)Yj9s!( zuGwOLVy=PUYTdGENmJb=H7g_O z`j|&rUS1b`RH;>A2M1;#afDCFqV+2htIcO2(@t%y>s@KLh1!`ka{cjxkNNmMt6+Uv z`LGskVdRuyg?H#|_^pGI-somXB`T#OW#N42C{f~@sgqS-h3!7oT!U?%`DH5MOoRqS)L>vK#>trURTk!!FVW*E*vR&y32yDGcC*>!Ig4!tXka+1M?o*t0o z96P(*u<CrKPppuX1(jv${}?e40a((9Yie?)c$K^%P3lJ(5l=GEVdvJqKD(5Onc%DT#WD zJ4RWdr~Xf-gbgqJ1F+kjJ-9`&5N}rpi zc!YXX)qvB&^y_4E_Y^VmWTCX+D4VnbHWaJ*(f}&ee(*NVMNDl|Quhrq_K%S`mJb^qeAs(xHQ}5z4a=<2>Z|f zGOg4uWRb;l%Aq@Wl#>JvrIDYR%Dnhc%WTHjEvO&7zEcHDd$-<)cv;R1GS`WrwAo(# zVc}zD@+ID^gXipW)BNdCjv8I6j9BYybgT%Sr5fQbs7F8(Op|tdMDQ*$F8&q4sH7=A z*yA!;cl)jSiyQLt9AV|Q9sO@L*xehSiU5L+fr`d3iM&DsAV05OK1v!BHQ84a>ox19$l=cPv*ZBe28*yI+zbSAZYLeOWAfL^K zcJ%_J4Rm6~dJXjz+W|{41G{WKA$p_%E%xd1Q1{+mI_xPF7+@A0epaUWfWABZ_FLaI z0`RcJ6ydE}99ukxOhWMPD0SGkWPyQlTMgq=qdFk;@KVO z#?S{de6qe0e6sFoA*lcTo|_l|^N26;R}?^c3Y(UN`4aqNK)O&M$g#(}Cqz%(4C(H` zl1G)!K|6yc2Zu`F?Ia5teK2JTlbH%iwoO!gg%L98FT1XTqcyRlDrC9Hd0BnAv?;xm zDfbd$e|^Tz!66xk(**?a3Ji6#S89v=~9T}{SE)n1gC-gA33Vj#ve@Y8nYem zl!`lD`dLhu9dfm&{A|^ComhRqtrlxMUgef1=glwa(DNJA`8U8@oT0}XLB|^%aE&e# zz?h==b<@RxOG##lWfk%RN0d8W$8Vhk6Z{vrhDxn_SxqdG_Beule56e*pR!W(1R^#>T3!lvwyONA3mEpi~4c+?g7;$jaa=}ueo>T)uO6cnpmlPPhPJyaqyEi5dNbvcDH#=B0)d@vayj5j)Gn0$hC7 z0ACn`j)>Jpp!BU(AJjcvHmE~6DIKq?>M0%krU+hBWmHT7%Y(8=-m)@U;WS<)l6^K! zsH}3Ene?4u3?1=tKj~AQ8Fgs%usqb4cAxmVIoD;@6eI_l3GU#dU3(A=DB{@dCC39$ zR6<^R&2^H-1>1~c8)BgcbOM=3iWn{*)9~g{5Wi`%o;QwG>&3{t7kctQR0*3`39q?y zVj=?9qK2(w@C*@RO<|I~kv%xX#xIJ*o2>G{%5T(PT>bg(K?ImNpn=affX>*L4>Y2; zb3y=H!$D0bN8rz|Hf-sEYWUf$dh$yd0)JWD=1z|NJCjmwW&u6Bvy(mQX)vNFAiWcm z>iSao`OAI_Ts3oSHx18&n!|F8^KHx5jbbx__Md0q)P&)B0!e0T*cmdXQks1Q0PZl{ zdGP%$)#c=&PO_u2(n9(TlbBs9Pvn(}R^Ma(Ef8Hdudcb3;G9KKF%E#iQ`A%G8J09$g)M1Jv+Rj~-awKc% z&(n4SShLoAb{ZO_tA0Mymm6dakeF-!bzd=coVXmYw|iU{o2`7Zr0+WVajyQ%`FsoL zy1q^NPT?%x-Dm>~CXrl-Xp3yH_rQvns{;-I-b(0gaxZP1PLQ8&+?!4@1j7#Q5ZMn* zp(i^f^rgQLpk{#q7V|!54bkwiMaNsGM% zycZ2h-uD{MS;ryQ@f2H9%>t7?v4E3X9g&t3Rsn(*RA+-G$adEIrQ8bJ|F{w?UM}U~ zDo}#kOO(Zr9+LO)8x2=PN;K1jG3&`*uMJa+1p6~&6*9J487cqd_=6&>wCOMwWXeX5 z{MzxG;R1ucM)~h)l!*!ymCjKjk3x^ek=94HiDoNWLz{Dmsb{dKHZ5W6TO+yOOzfU; z(`3>_$}7G%Ds005LD}x-jaFdGa=@Z8ZrOf8{(PW8PHMdYq3H!(@=08?Mxg<+m1+jM z-+adP-WceDc<&f)=klU;5zDKw4WP%Ato7tQJ{*3TOu5!znD(6cXl}A<2uwQj^>Cf! z)SNF-E1v;fyhS25>WSRA>c(XNGk|q*r2-rjGAYRgwNwhI7;|0Nvz5JXozJwz= zEfMXtRCFrk!$s7L3e$2vtRLFrX3*MgY^|YZf&2dQo0^v;`GxH#>Ynb`f#ND*_w}Fo z?`e1tvG+hhX708$GF>z{w>#z{21-u5N5~cjnQ`%>_VNzyt^s}JZkN)=FWXB4p&M?B zp!}YliU~hrF-`?)vC^sRv8s8!&&{KUMYNmxZ_!VBIAs?yk;CnjjFAuN%`0iuPQP#< zSITEtl0vegMAZOW%1GkX=&y1!Cj&^`4alS}nEsNHfU&l=wt4>uBSdXtbwyowZ<4*I zT+Eum#F90kv_T@%8|$M6Epu$!*dR3VHKyg|YdfVcoo z$-@Ai&J&5JpvY`~;qX}5-kY#x1o$-bHnHts8&uoDqM6lt{J)J3y=C7&Y-_uYor^es zcObv7gbQ@7T1P2{)5L20y9PJ&xF8c?ntK3??o_)^rrxqWKUic!J>qV3x>C`*y16H6 zc)mwjS(rmUNcMPMrp#V@ge+ez4kQywrOKiCOro$&fj)+F@jhte#tpQrQzvvJB;51N1eT3#i zY(IRfzE6E57D_X)VOiJGB^Eet4&R!?BEgLA$D3W^%)tWQRjTKQi zy(MNy^on8t>Z&H?5*Y2@T9Q2s2kdzMWD#fuRT;(lpga$c9_lX)* zwo-6Euxsy8C|PbP=n5$Ur*t#=m%v1XMZQd70v zw@VDL#cNgWlUivK_9NSGlj}%zT{Ze5{HEnhS9dq5znKVFE+)HpMIsX5I_v&XPd1P7 zYYtiQ_>rdJ2LBQA)$IX%<0`bN-^p(Qhm1O`=#8$w_V=qFYu_aSE&na?3c_|M`mryr zicYj9tb;qyH6}X)_SBi+RS5OVkdBI5%Z=cQwPAP0qSEsL7|pbCaPbs(@cSv$y)%nG(*eX!G@DS+2v` z>5TJyu{mOU5jhuxf%a9pT<$J4O+wn7+=J13FSNN;-+il&fKO+*)7o2Axu&{Zkrw`- z-B+4zpOpgcHoT-&!NHyaiN`xsLaTWxKjV@X5vx z5YKx7nDMQ>&-ymIU+lQ*}>>?TLiFc z&;!hgeCU7wPzW&36eRS-6h;rm$Z=xP9rbJOX9qyJ+@A$X8WMvtHkOdneg z?5%ohCuUs9Bwxq0wpb$(8EcmvjP5Gr@{TNV94s`9dD5LeM^x|uoJ1h=-MZrbbi_o^ zR(H_hQh}YFt2`dyu?+O#`_FHIHegeC4MUqZ`I>yf&(<`_eWXlWB2!9GXQv0pVVo~y zA{tK^ucS*f>M25tYUdeg09@-v!|r^kN* z*YQAu3nTN^0G%1|r*NA%j-_Gd2edfwnrk@>fV_GFXBx7OtH_)kh==Q@wGu4dLZy*R zvzv1>ZFOsJ*DCZcS5`Z-Rt>uj3CBY!ZLNct*vQfq{KwnwC|fza^PCmu?XUG7S^AJ7 z^1yZ9B<~$yT`xe_z{kaV>u9Coui&b$e1fvdVdPdw?ZgNzxcxFnAy(+i>IaTFA)HSq zees>PvtOcA9@R?3&&SZyS!ZLYz*yr3Z1Zc`x^80tpEpahyQQhs4k&q#jy#&KJW&q$ z;N#{vz^_wK?ZH#xeSZTKs-)%73W3JJz4>rz?rTCqB#k^42Ap4Oj*Gq>CxUH(yO#+;ujpiem6d05oeN*$vyS%=6Cm=2$@#X!|u zjG7UuHGHn$E`h5WY{C~phZvPF((SvBl+6x5EIO~j$;9R__v?b};i6wi(hq$A#J__%%m z$8C51tap9+Gg#*QFx3*s^g+kRTgb*1Kuxe80n<8A4g@o%$?fy>ZLXh7Ls*rDJ}gP+ z@;#3d#gHim?&fM^MZe?Zu@%i8^_O0Qs%X-8*7s)+fjVlX&RJ&{jZ}+kvt8;J$DM6- z%#;^NNUA2C+3bx?eq6}9Uu2wm_no1D!k?iiRR)?~|Bduxrk7e_fI2-Wf2=m&o84>X zfWx3NKV!srL{rCq+-uL*Mc6fwp?x@n^7Fad%2NsOLu&xUc~B4%F13-BAa4$4$1*6dI~J zN%pwB0mRkdt@-tlk9-Oe5-KzW28eSwOr_f+t4@zoQ+U6-R$g2BynPox?jD~RQ(b}~ znEpcl>!RQ5d+dkI`U-dRv@j!pUi$O`#*XvukrT)<0uMG^vRot$c>vPd1JT33HyGUL ziNAd787S4)TJ69Snc#`Qm2zqa+*7th#a&Z^H%Z|;T5DEz8~X$oO$DRZ!h7tnJ8Oi- zE0Rwx0>e|;Ot*w*0=Y2PqV+D?C_ZKXs4qjxNd+`=e!rdz$gbARrpo(@@~X#c!d`)& z=Y0l0FSTkU383KTd4T2p1f_9F6TCg;lJXn_a>&eKgT!Ef9OHj06H;jUu`isK6^h*a zCP2~@_Esm86L{hMdKG=m6{sbHfZY|iv9WRg0|g`eC3Q z1HHg@5X{v&l8n9zZ(9O-R>;Mz4fhAL&cC8QP|e?2)=O7U`ed|=Vy*#}0h>>?l_lMJ z1SMm#-d7)%a1V%me`Q_JH2jX{a4#Fly9B}L6%4>*Zk7O6INIPZDP;U%8pyb|ga&`O zh5=2reHR;vlLM)wDo{yMGC6Rq5exJk{uCX9yA^#~tw8wvs>W`ag1WfmKm^_BUoENQ zrO_RQgAB$xA7bZpYXldVT>rkVaA~fQZqs0P>tiwE5WU zaF?k4pNo@)ujN#*H*WcF{%(+{{+fsG0RQv+#)mC)9)yx`JBr=PF-e;*-f)$4gZJpa z?mPZ(uK1dSA8dAs6w1-Z?_v!-eieii7yz)}vKDS|&s zt+RnHIRID!&_b~&DR^vbR0)MI*nzmuifOuxM-R}4fc|^c%*+hlfzEvg21H0O%Pm0Q zqk}*)>SQb7@Qlg|Aaz-|dUJ|eNMxs1=_r0;!h%iADig&Ck45vMSl7WF zkkZ@B-nA>AMRTaVyj^>a4863356n7Q3QW*cstGt;%Gf%{NYE%8b>@QV$0x_bAy7ok zjvyS0xMiKIiMU2WVS&iMY=T8_sdlKOaPggEXySs>R-m|{G0!7~q59|ld+GmE4Ya+{ z%2G%h7|@!obxna^gX%}7=(AcyI6L8GgCefcs$rWu zevcOioE@z;{;BH6PF8a2D&yfJ@m$dON|ayt8UDlv*LmsR@29Ebk78qy45oOlWnPL9 zM*rVMg)-|U)extNGPYu)(x(HdcgPA6T-{$?u@~ItTLQ+4jpIbx4qxy;{sMt@(;HW~ zO^F@vt*5}`p!$LuTx>UdbrtunrAxU#MTZ;=1_@>A7*yFolSQnC*n|GWTN78Y<}&7J zr%P3$Lr=W$dlRE|?lYxI6Sc0!fxFmFdSmnWM_`qf;O;pI*uq^M`3;G-({0VZ6p{Av zRgbnKhbS73p|{t?uWIl@V_rvEK9_s;`+|4(1k%3%6+aR*A*pC!3K7CT2rqV5de%BrAAsvU<{a$ZvZYFLb2 z{7R>r|D@KXM-HVRB{1**V+2E=<_(`GnqFPo@sSaMlESWc}mTXA2givWOn`m;J(;a`9Y64PJeU9Dk~5PS@}aj6I$!1*6z0=G z5l@#AgGgLa*8M|XbIsBCd2VQ*XDkR5N!JC`edP2myW+xDm@e^4QVhq!+rNu#1SSpn zXOBXc;-S36*x);`SemrO%=+BSjCv#df{F#9!$ow?HVXV+`J+a$flP{nxfD1vaO(_j zbaz$|^^gpw8-=_N&sb3c{1#N7k~MehX<_C?Hn|9Q8{t=1 zQJj($JJ!tQd?Kr&Yx4Rks@T{hMads#QYZ`^*#~M258PY$X@=ZkSF~OmvD;Bc#IP66 zQoMf0f~i0*Wcn>68Y~>GCICh7Yk-Kq%B>kybTFXBBnq+gzxhvwuPP!Bn580=&+aH} z&9k)~^G_*&V<7nsR2=YuP9x_XhNBF0BhT>v_8GSDnu(YFT7rZu z5WAY-_$r$<_zXAQ(A+ZvU(F*69NOT=;&UGeT~eKpx~{eyJ0G&;hlbYa~1Ceo^h z3yQA38uMykKy=UY6Xjbria4bYM8n%W&1JlKiS7IWnj;GhBW)psRt0cDBclyARkeUI zT^T>krb~qPCLrJW{7Vsmv`~ zO7v=kT`r^EnnXrupiDe~4;6p4rmCPIm2-FV9o*tidXc4QlO3{vq`@?;6LqP9;M~3C zvFjwWKZy&T6m>!?5!GP_w*EC;6JA_KJ}lAz`7a}Zr*R$AwKilkvyHGWG7!D&3%8)% z?>KiaC-RYooJ8$mUta7?&?N?PO|dC9-U$WSY?W4MKXyr#7U>01Z!ySs@M>vAWu{Kl zfl|VQ^=F}@J$I%CqDBC0t^V+#Oox4(qV=>y1Ad=)7GHWbEN$RlRhgqDnC!v#|8$OX&5-U2toHE}g;Jii%3R*00K zJ&Vkhv%PZS*5rst_R5Fv7TQj_cF+X&l>w+!99PkT%rbl9>}*B*$@ln29hWY<6W}hpnb~zX)dT&CF;%VI(JMIOV;S4(rl|(nSS999TH%~=5oB~-7C7MY!?b~S^sDSl`i$BXTK-;~^)(}= zs-tf3c(;oJ*2jlA5-7jh>`M&N4Fxr9H>#KS{=p{ZXvz`U75Sg1V@=O6a>oN{uz@?6 z@8egw5jr6L^X-SG3H2iwd3@NC@elO7FC->#gd~+ z;lHEpGN%~k0Ix->;1<>%Di;;t}4`AagqC(1_(bxqQ#MeQ!b z|D?LHLq~oJ!tpM~<|KqewaHui8}*o+m;b~Fk1VrtM_zQAllyaES*SA^{K>pA-vEl} zaK%njE4#^64=aiK31yHep1edFfvO~`gNoTr30PDWE*Vu8Ba}#`imNtQHnc6&`r;>8zIQEBFlpP zI!J8HCQC|S$;U$VMC<-?%q_kv;mS}KeWBsN5<4dklFl`G3X2OXf2!+2h5g%&2?Un4 z_BY5NH}nwB2y~IjE!8im$S8Ff1Zvg~7KGC`NV)=7N|mE`8Gj)?e|cFL+NZ43SlqEu zXD0LX83B($(jcGg?X4#OC4tDS)if+6M@B{F8CeAh{oz@ieb2SbL1EFX)fy;j4Y70R z5+Cj$;E9IeV`pgcEvj?RbsWZ>4$0>lmxJ=6FDMA#@Mg4Xv6r3UPU-lp2V6Rz$v>s} zYP=58)bU<>=%X@olRV4z(2*h#QdOv=)V=jw3wb7G_WdbdLC2Xe0#WcEcDVw2=unor z*O~H+zsJZp_hRUY0;emMkD1hn+i|-Y#u5*8|dmk-x`$q zB83YtmQuQrLu?(6*9u8tD`}1+^`7Vmx^C{wC^Vj$Jxsyem_v#6N&P>OEE$fezgV>^ zZMB^UL+tw!ZPrE#5(SKfF@|1KEwWmj?YmbKg(z&lbzq`-fiIEs7S8&09R6M0>0`Ct zn66bg#^#an#27d@OfuI!b2VvYp01CU^&xtBWCG_d8%Y01kUi@HOP_Yf1OtI$M<_bJ zJVt$@(n`Bp_bLSc?6ajJhufytZqOrkfe7g$_AO^VPR5pJMrYB zDyWG0UEKhByNrLK2c5&r?^;8CRk_Sh#GUu29EqJlqOivshRz>q$kQVm^Z^;2^{^t* zRFN+M`rh52@g%JPt;}+(FE%&&oJ}IgFWFp(Vwp)t4kpFG{>~^gDrVPZRc|&Rwv?dXB^djPpD*efr!*QchYwV-j20Y-%iL7c^egw0`W!@Tl zN91HBfYlYDA;1;+WEq-bN;qRL%YP>Zh5>%e0OQz<`!tH%o4~CVV)E$p5x4Afi`C_i zotPAu09xc2R@gbwC+wRY&Y*LQN@1O?fqS@0qf*>)J(?%O9I*{D(XH}luKq{M6{{x{ zA{}`5hlqJkNQ?YN_*$%c(}! zl{23bFtxR=NTKcIrhIj<=wZ(aTI{o(v3`ENwzKR^mg00H#BVuI0A3qzCxr{wGbUB$% z(TF^FTw)^H$mPsfiK34-BoEx3tURy`Sl%bxx&e)R3q%6BeR5K9aRP7R$@5a%v|u5_ zhg>F=iyq|qzn>r7bYo1>zz^%AwJtjA!%YaHAyF-@uG>(3k}4E^%3|zG8)O#-TQoLY zp9(EdFVuJYK9H{%!qi=m;WOmv|Le6LG(Ran*fv){3Am9}FbR$Fz(n&_+uka*!uv5K zje+;t(<-?Zvfy)LzYv<#n60q-F>ZPbaB{=kbLIiu^_ih^Nzc)c>Yli08s!Lq^85z{ z8kDSu#ySdrZ!Re2FdJKpzdG4u$+UKYPJ4mZk~JzR4=E?jbGwZEs#9+3y&xfF|5Y=3 zM$$s{r}v#S`;bumTWy|Vu8s2)7Tr%A9C+}@_{{3zP-Z65#Kq0#bg$(kx~n=o(3u_; z3FW^k2ww{-GK>j3UEDGcSU%}+hfs3Nzh4scp5|&QBQy&94>8dV z`_qyg-5m}*x=Inpd&)a^rJ#`?B8Bk}_>+E}oouZ}vG9}do^$ozY9*CfdM{$yBC`fZ z)VMP=)lh*?1Z}mZqnf$e<)G7DDk-6n zW*&`^;4{=tqzd;vJa_;|ny|{|IdJ6dP#-CC%++ik;Q+8=vz_3v;I`lk*=$$thVmAe z*7E;w0w7*;`MP<^HbC2(Xv~IQay)FzFy(>6IW?9>@3wyv4Z4L*9 z`jES2M<)hJt+^A`+^wPUDpGjkxj z+E3=TfUohY7qme>kcw2<+>ZgUm0c0mYyxdg1Up2$5pH@_m=zdf=t2b}hY3L9l*rv@ zx{R&3nRWH=09(kfo$d@5`fbhI+2uc1CgFWlwIuf~F8bR$o8My*q#FDTf^a^LjF9~) zg`EL$J3Bkb%c+J+THweL!^+VsgMXWv^3MTH94);SA6l4l%;pvsr}K~7YLidJ@27t| z-h*KxVQU2y9u^~f+}(XDDce7a3bQl_XLvrnkLPi0R>59O(d@?xMTQ@GIB4?rbd{iT z1^G;-Sbe4S=~@Ai7Th^*&E)eV10S9-52mB<7ic(t(x#XgkEYfJ>n2UG{8x?2%%Sxj zbGOav1c9xjiRzY>l{F998;7dRe4^8l2~V*Md|gr^2;tFbiR951E^r<-PZo*R$y|tH zF`wq)BM^5fYg}RG(f;stq4QzQkhiQy=Yuw$ zktsuS&&ZFn!o@h{ ztj?F044UB1kSOYq*TYKdD(#^xx2@QDbfYd(S>qP$d^GrvLFEqwG?OE!XA>Nlq^zj=Ue5)1GeR;wS&z=)$k2urdvOPswPL4 zsgyr_$DAl=im0Ct6xgZ-7$BBL`Bb#iaHd@K4MWRzp8zPTq`2S3Kwjq60LA38*O&FN?ez!&2(bJ-!q1M{-x zO}TIHSb8#KeHs>CP<;RX{c!5DqB4HNGb}G~@Asd&84vi{NF=?0(VyMyfZZQz;z07g zB!NpEfD&^?w<)S|w8r?MsZQF(f1&UaN70Boa+{yk0JD z`;C%Vtb>oL6OM*=x6s+Pv9eLX8eVG!DPj^FelKNDamee6ffDrmh7jwKd?>5m ztXGM7;|tK=UuqgFr?Cp7vS{i!TD`BNbw#~<{>$TH>`c9gr0dUoU+4{u@D)a3bQW}N3fHQU)`>Z$bnj~PnvR>XkG_f_B52gDS%cj{K zLr>zr_)gjhcpv~wLVLeE6Zv5Hz4scRNJRfQE-qhjJzr$PBnQ7eH8ll%!~od&A`QD1 z2=jJjEk1THy!*A>UHT zeCN~m9RuL%*X4cP=Fy9#_46MA)eB(YTm6Cs&{Vyq%XQEO#4hX|-mm@TR$KqUZT*@j zb#A=_Y2sxB14?ZgWHVRmNb;Ijz}W9#*KKmiXocrAW5Zxn+7x$~Z_{|w$4l!i#SUKDri zqBo8xCrq?=IfZ&~N<+5!GZ(E$SddvC7|h{ymdgE~Ql-776;a$`lLR-t5$86<3S9fjEp{tt@2 zU%&J8JTZ9_f6*MFGM?>tDbkCSH+Xq@LFxcyBQ*bUu~CUxeJuD^92klu25mM*gM>N? zIDd=eTp*lQs%SW91xRoJIV6BKqbri4{OW4d1{#j$PM$E;oe>hh;8)mb&aGKW}cCL^Wh9aUBW(~l?5DT1X;60Reppic0`xh zCkUEUJO|Mz3dcTj8T%Zs9t#SDIB_op++B#$F_!a5jZ@`dffPX04PK~bze&2#%~{yb zzziHIh2(=)mn6wjv^R@p5Q>!a{9;BfKO(X)t?a3usec7 zFJWLj2BY2r(V%=620L2`+`Mq{ygr2su0Kfzi!jf1q+m%!4Z@ zCLZt3f^_6IU9)yVaG~rlwhHM68SK6JXBB|jV%~!Zm#D{VLyNm7qll%}8wR0NdB07O zHqV1O%zO~0^h*^d`@D+-)34}F|e`OocoZT+j(lJgaM$7HC0hrL7TmxE;-Hm_b1x+Ckq9vX8MJWpTx=WXZ+gV z?#BopmbXW_mPbCRU3krFz4+~gmWX?HBt>qJpVNGc7>pseDT_1l2oI_8L@PBukBo>E zrq=^S(ry>~1HgrA&&%7-K7ddFmG7<9-kBkJ6FhvX8dFr5^Q0E`V4XW9;@L)Ppy?lt zp9pTu>>SvsGe{zm*on!%X}wHS%q160nN>|9$?2ISU9PDC_kaL_ z^FS18WybX;JJXkH-fOHaCrq=KH&_xnwTX5ZFGiu>`#|fV z*_=l5$uyG%kGho$4fz_j`yQburuhNZQV@6SKnr z`BNB^B)^;e>_^qZCXi+CIH7Jh3}+E@@V&j?<8*@-CMH<@^Spz z9g$Ox(T^G96d8Oa-R2dp9kFA49jjw~xbi-r~9jaP8F3}sX$Cu`tsg&j(_ zdDgz5KPLVgl5c_mrK|s-FIlJ`m3!ZL^b-}N*tmS?WER9{CyQrGwOs zV@^&-26NtFAoU&qF(gFvI=`9?axL3gb#>K3%T8Ba<^AZcJ3Sqrv@wv>PufWRhSjW{ zD@^V)A5~!|8x@2z%$3|Uaiqwbxn^WbfK6D$ma&b?C8F3wyo3FbcLeB@x&S(JbX!P0 z{~=p^JSYPZ>4^Z<6J$~Zua`Z`@_Kl9=_J|Ay7NWmu~f`mIkr#3y7E-iGLX#{x@6(V zkbN>1=8!Yl(?t%aia-nQ+}%b&gsIheDNYqM2jl0n&fb`Kg`by^XarG#AO zn6tu^2390IsqrV#Ct@B=WM}01k9h&P1ooxT7auy%ICKI-*t4G(+c;U>fJDh6~ ze>Dq6vtnB~0aEPt0z5;9d;Q!uHow_;*eyKJ-vw6xg>W$ZPXpbL)alYk50KBVOE^7J zS@Uh3oi*pt$@p&JH+P{U=iW=uKDmX5p#yfhsNuTGY!BU(okqAaAb!)%Hl@umWpL(C zx-@>Dl@H8KF8qDTquup76cq5~KJJuC8E`e>eWYcB1+4SJ?(cEe^z3vodyRSDB7T1s zsfcuJ;6W4yi5wHS>%ddvOn26){P!#DRz^(J0h(qMJFLW#Mac2CmsuB-A)HoR({4_l zK+zYZd$!4X;`iceJ<;=G=x)G`=O%;uyoXdnva1={ z??yn**^jtN)>YlZi_!b|VDeYB&+4P?&EIKnwbK^{rEPEv zFx$mw4v!=AUIeO6h)8NL^h_mnbpoCbjG{3I09YFo8wrIT(t6If$kOrZXvrjn@@Gah z`e3>lk&APF^=H9zO|a9ys4|NN%_hSlJ8yedO}>nQ;Av z2UR`ny~|^DU_Mg=FvC626z=G-H!mlfUDc^s$FExvegQXM0uY3k)3%?!ZYWc}x`#ao zfQPN~qk!JHy#eM*bu(SD@cU5_@tCi#ZNB&oWmJe;OjX(Jf)2qSDC%doYfp(eSprYw z)v&$sY=YS1%BP2_Uo=y2KhkH&r70}*SOBVD%?c6*CE&sQTMaMoKE6WOa+0x_$E?hg zdJn3IS}dnptQ{!sWQ8+zA0&yLW!!TA{U0sc8ccVDO|N>CuqY`6Z?FeN=mqVpVEhNa zy+f+&&+IR!*pTKx6YAV3cSl2_0P!h)3uJ~L z8Qg0DsI*@IjF2xrs(jkrLJ%IMaY4vG9GUF$+WtatiB(k5uU$N*`1;@NiHK{@nG;Ir?&WTR`PX3;0Ch7ur4ErFis z8?h_1jZ=IGkRY*%D$EW&mNF)}=1^pq`n23i`R%;42WX|bYEHkt#vQK%rZTr-KG>)K zq^cSj#4PbVaIY)RW6r;f#B-iCwmxL5wPN^8ek-Eq631!qE)@ornHNhQ^taxLw(j?1 zRzRdV9$V5iXj9lJ^t0Z{#z`~;nkp6Iq<&~LNCQTHTbaVEuv z{95G`U+ATz!z*@O=37J&*^(3v!_;vXD7yBRMA|7iPqO(<#aM4tZM9)+m;21v=>=4d zwjbrHn3O%P6Z`}CYS|LQBG{m;xX=2c^61)od;GCbT#^7%ghjx(D~8_2*Ji2;=e0Hh zYwrM{B(bVMGwf(9D~WzzyWC8A;Nj_hJ!Vnq#GlgH3m~ZOyMXAyjup}QtPEgEIHt>( zZtFemXWLb-^R{Q23&q{LeR{Nm0RuT(OqXaT0Tk`ZJTdY6@vAGFvuMnv=8cQvUm~wZ zFiUpjIcR^jkZ)NTblq0F%YPBv6jB8~1`-RqV^)B}h}Ts>RW#U?>E&gIc_5=i#Vi{{ z6&OUQgMJWjhx;=)lge{oKr$O#7x?V<7RO0}fs42}vKH^X++}+(hwQbb7l?{`{N?WW z7c=&gEe@qsC$sji0Gft_>tF!0C&J&}x+x%V4ZA-pn5+}+y!(y?X*>D}_x<~2?hK7c zT6mGZ*{9`#)0Ma+u8n*{{y)x^P4E$@VFvxj6B%ykNNm;rzU3G5ESIpVOvnXYKIo%E zizfg((Fd|y!XthtqRy&Pc-KHjHgYgsTHGJE&@$0>cI+BO&DM9cnw@P8E-2OMf)!%+ zegEb{08F$@9Pe-Rhn{Sx>FDTqcXC0eb-J5SB*JZnPl+`KG*u+n@nS}RTiN;iVDGu; zZ+5q#Sft$4L1-Qghn*pJ*GAJMT^Yav%aKa(y0C{JsX6?Issd=W1K1G5E`DU7hp*^)J2V%>qo2u)16<7j$M&{J*V%Ta{q`?9sHw1M=KGqY{&NI4iU{HL~rcM&j+t z!Bu^IaOLpMX6u1vS?m6jF970ZpY8I#vEY9wxXT6~(^Jy}90R{miF%5#1tOO7M~SI; zRaF(lq%Bw$Fqg&M!t|U);HssL+ylL&^z?qfG@1erAr>*8)gf&;L`);`7NhsF#?KPd z>hX1e%Zc6znkY_i50LyrWEKs3*r}Ir>+VRtt$Mo?0e`RxZUbn_MdhFS)r}5D6&dha z`eQDrVKrt&-k}J%Xl1TMjghj2%N8+T^Pve&U?2%X;(L&LoZ|k%2iLM(wK(=f5i6xm zIM|&N&UY1H#6AI^H|G4IAopHPu?m$HFe=FbA#CoCiDYOr9FPICaQ*Y<(f>#1Bp(3X zFOhI;1jpO91jvir)YSAXvm}3~9pW0t1Q^m&aVTX7TQrFVY}>g^lxYKAZ4H%63W#0U z$PX?)elUZcT~4tY0?=IcPX*E+M6UnIMm-;hIuOAB^Y4I@D*`QDcG#J$L7~oaeLLPA zTuSAWv%GH{o~|6G1e!Haz$SJYxQm|4s2GP&SJ;pO6&Kuld2y`FA|1rud%$q&fuPOL zgX$XsmizO1S#vC-rGL2Jd_Xkll>mDOZS-CS6r$5VsfgDA0MU^J$e>iEkmCwK`y|B2 z-^*4C@dr0g*QTl+NB#ogV*e#V&i}p?{r`mqa0EKmR+F>9qbuOvJOrw$^Po({`o;eM DA~>$C diff --git a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb index 4d0f96408..f34180b4c 100644 --- a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb +++ b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb @@ -100,9 +100,9 @@ "outputs": [], "source": [ "import numpy as np\n", - "from scipy.interpolate import interp1d\n", "\n", "import openmdao.api as om\n", + "from openmdao.components.interp_util.interp import InterpND\n", "\n", "import dymos as dm\n", "from dymos.models.atmosphere.atmos_1976 import USatm1976Data\n", @@ -158,7 +158,7 @@ "## The cannonball ODE component\n", "\n", "This component computes the state rates and the kinetic energy of the cannonball.\n", - "By calling the `declare_coloring` method wrt all inputs and using method `'cs'`, we're telling OpenMDAO to automatically determine the sparsity pattern of the outputs with respect to the inputs, **and** to automatically compute those outputs using complex-step approximation." + "By calling the `declare_coloring` method wrt all inputs and using method `'fd'`, we're telling OpenMDAO to automatically determine the sparsity pattern of the outputs with respect to the inputs, **and** to automatically compute those outputs using a finite-difference approximation." ] }, { @@ -198,17 +198,17 @@ " self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r'])\n", " self.add_output('ke', shape=nn, units='J')\n", "\n", - " # Ask OpenMDAO to compute the partial derivatives using complex-step\n", + " # Ask OpenMDAO to compute the partial derivatives using finite-difference\n", " # with a partial coloring algorithm for improved performance, and use\n", " # a graph coloring algorithm to automatically detect the sparsity pattern.\n", - " self.declare_coloring(wrt='*', method='cs')\n", + " self.declare_coloring(wrt='*', method='fd')\n", "\n", " alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]\n", " rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]\n", - " self.rho_interp = interp1d(np.array(alt_data, dtype=complex),\n", - " np.array(rho_data, dtype=complex),\n", - " kind='linear')\n", - "\n", + " self.rho_interp = InterpND(points=np.array(alt_data),\n", + " values=np.array(rho_data),\n", + " method='slinear')\n", + " \n", " def compute(self, inputs, outputs):\n", "\n", " gam = inputs['gam']\n", @@ -220,11 +220,7 @@ "\n", " GRAVITY = 9.80665 # m/s**2\n", "\n", - " # handle complex-step gracefully from the interpolant\n", - " if np.iscomplexobj(h):\n", - " rho = self.rho_interp(inputs['h'])\n", - " else:\n", - " rho = self.rho_interp(inputs['h']).real\n", + " rho = self.rho_interp.interpolate(h)\n", "\n", " q = 0.5*rho*inputs['v']**2\n", " qS = q * S\n", @@ -457,7 +453,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.0" + "version": "3.12.6" } }, "nbformat": 4, diff --git a/dymos/examples/balanced_field/balanced_field_length.py b/dymos/examples/balanced_field/balanced_field_length.py new file mode 100644 index 000000000..a6dd223a8 --- /dev/null +++ b/dymos/examples/balanced_field/balanced_field_length.py @@ -0,0 +1,189 @@ +import openmdao.api as om +import dymos as dm + + +def make_balanced_field_length_problem(ode_class, tx): + """ + Create a balanced field length problem and optionally set default values into it. + + Parameters + ---------- + ode_class : System class + The Dymos ODE System class. + tx_class : Transcription + Transcription to use. + + Returns + ------- + _type_ + _description_ + """ + p = om.Problem() + + p.driver = om.pyOptSparseDriver() + p.driver.declare_coloring() + + # Use IPOPT if available, with fallback to SLSQP + p.driver.options['optimizer'] = 'IPOPT' + p.driver.options['print_results'] = True + + p.driver.opt_settings['print_level'] = 0 + p.driver.opt_settings['mu_strategy'] = 'adaptive' + + p.driver.opt_settings['bound_mult_init_method'] = 'mu-based' + p.driver.opt_settings['mu_init'] = 0.01 + p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' + + # First Phase: Brake release to V1 - both engines operable + br_to_v1 = dm.Phase(ode_class=ode_class, transcription=tx, + ode_init_kwargs={'mode': 'runway'}) + br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0) + br_to_v1.add_state('r', fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0) + br_to_v1.add_state('v', fix_initial=True, lower=0, ref=100.0, defect_ref=100.0) + br_to_v1.add_parameter('alpha', val=0.0, opt=False, units='deg') + br_to_v1.add_timeseries_output('*') + + # Second Phase: Rejected takeoff at V1 - no engines operable + rto = dm.Phase(ode_class=ode_class, transcription=tx, + ode_init_kwargs={'mode': 'runway'}) + rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) + rto.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) + rto.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) + rto.add_parameter('alpha', val=0.0, opt=False, units='deg') + rto.add_timeseries_output('*') + + # Third Phase: V1 to Vr - single engine operable + v1_to_vr = dm.Phase(ode_class=ode_class, transcription=tx, + ode_init_kwargs={'mode': 'runway'}) + v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) + v1_to_vr.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) + v1_to_vr.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) + v1_to_vr.add_parameter('alpha', val=0.0, opt=False, units='deg') + v1_to_vr.add_timeseries_output('*') + + # Fourth Phase: Rotate - single engine operable + rotate = dm.Phase(ode_class=ode_class, transcription=tx, + ode_init_kwargs={'mode': 'runway'}) + rotate.set_time_options(fix_initial=False, duration_bounds=(1.0, 5), duration_ref=1.0) + rotate.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) + rotate.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) + rotate.add_control('alpha', order=1, opt=True, units='deg', lower=0, upper=10, ref=10, + val=[0, 10], control_type='polynomial') + rotate.add_timeseries_output('*') + + # Fifth Phase: Climb to target speed and altitude at end of runway. + climb = dm.Phase(ode_class=ode_class, transcription=tx, + ode_init_kwargs={'mode': 'climb'}) + climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0) + climb.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) + climb.add_state('h', fix_initial=True, lower=0, ref=1.0, defect_ref=1.0) + climb.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) + climb.add_state('gam', fix_initial=True, lower=0, ref=0.05, defect_ref=0.05) + climb.add_control('alpha', opt=True, units='deg', lower=-10, upper=15, ref=10) + climb.add_timeseries_output('*') + + # Instantiate the trajectory and add phases + traj = dm.Trajectory() + p.model.add_subsystem('traj', traj) + traj.add_phase('br_to_v1', br_to_v1) + traj.add_phase('rto', rto) + traj.add_phase('v1_to_vr', v1_to_vr) + traj.add_phase('rotate', rotate) + traj.add_phase('climb', climb) + + all_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb'] + groundroll_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate'] + + # Add parameters common to multiple phases to the trajectory + traj.add_parameter('m', val=174200., opt=False, units='lbm', + desc='aircraft mass', + targets={phase: ['m'] for phase in all_phases}) + + # Handle parameters which change from phase to phase. + traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True, + desc='nominal aircraft thrust', + targets={'br_to_v1': ['T']}) + + traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_target=True, + desc='thrust under a single engine', + targets={'v1_to_vr': ['T'], 'rotate': ['T'], 'climb': ['T']}) + + traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_target=True, + desc='thrust when engines are shut down for rejected takeoff', + targets={'rto': ['T']}) + + traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_target=True, + desc='nominal runway friction coefficient', + targets={'br_to_v1': ['mu_r'], 'v1_to_vr': ['mu_r'], 'rotate': ['mu_r']}) + + traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_target=True, + desc='runway friction coefficient under braking', + targets={'rto': ['mu_r']}) + + traj.add_parameter('h_runway', val=0., opt=False, units='ft', + desc='runway altitude', + targets={phase: ['h'] for phase in groundroll_phases}) + + # Here we're omitting some constants that are common throughout all phases for the sake of brevity. + # Their correct defaults are specified in add_input calls to `wrap_ode_func`. + + # Standard "end of first phase to beginning of second phase" linkages + # Alpha changes from being a parameter in v1_to_vr to a polynomial control + # in rotate, to a dynamic control in `climb`. + traj.link_phases(['br_to_v1', 'v1_to_vr'], vars=['time', 'r', 'v']) + traj.link_phases(['v1_to_vr', 'rotate'], vars=['time', 'r', 'v', 'alpha']) + traj.link_phases(['rotate', 'climb'], vars=['time', 'r', 'v', 'alpha']) + traj.link_phases(['br_to_v1', 'rto'], vars=['time', 'r', 'v']) + + # Less common "final value of r must match at ends of two phases". + traj.add_linkage_constraint(phase_a='rto', var_a='r', loc_a='final', + phase_b='climb', var_b='r', loc_b='final', + ref=1000) + + # Define the constraints and objective for the optimal control problem + v1_to_vr.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.2, ref=100) + + rto.add_boundary_constraint('v', loc='final', equals=0., ref=100, linear=True) + + rotate.add_boundary_constraint('F_r', loc='final', equals=0, ref=100000) + + climb.add_boundary_constraint('h', loc='final', equals=35, ref=35, units='ft', linear=True) + climb.add_boundary_constraint('gam', loc='final', equals=5, ref=5, units='deg', linear=True) + climb.add_path_constraint('gam', lower=0, upper=5, ref=5, units='deg') + climb.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.25, ref=1.25) + + rto.add_objective('r', loc='final', ref=1000.0) + + # + # Setup the problem and set the initial guess + # + p.setup(check=True) + + br_to_v1.set_time_val(initial=0.0, duration=35.0) + br_to_v1.set_state_val('r', [0, 2500.0]) + br_to_v1.set_state_val('v', [0.0001, 100.0]) + br_to_v1.set_parameter_val('alpha', 0.0, units='deg') + + v1_to_vr.set_time_val(initial=35.0, duration=35.0) + v1_to_vr.set_state_val('r', [2500, 300.0]) + v1_to_vr.set_state_val('v', [100, 110.0]) + v1_to_vr.set_parameter_val('alpha', 0.0, units='deg') + + rto.set_time_val(initial=35.0, duration=1.0) + rto.set_state_val('r', [2500, 5000.0]) + rto.set_state_val('v', [110, 0.0001]) + rto.set_parameter_val('alpha', 0.0, units='deg') + + rotate.set_time_val(initial=35.0, duration=5.0) + rotate.set_state_val('r', [1750, 1800.0]) + rotate.set_state_val('v', [80, 85.0]) + rotate.set_control_val('alpha', 0.0, units='deg') + + climb.set_time_val(initial=30.0, duration=20.0) + climb.set_state_val('r', [5000, 5500.0], units='ft') + climb.set_state_val('v', [160, 170.0], units='kn') + climb.set_state_val('h', [0.0, 35.0], units='ft') + climb.set_state_val('gam', [0.0, 5.0], units='deg') + climb.set_control_val('alpha', 5.0, units='deg') + + return p diff --git a/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py b/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py index 13b2edd6e..14c2b3caf 100644 --- a/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py +++ b/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py @@ -2,11 +2,14 @@ import openmdao.api as om import openmdao.func_api as omf +from dymos.utils.misc import om_version from openmdao.utils.assert_utils import assert_near_equal from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse import dymos as dm import numpy as np +from dymos.examples.balanced_field.balanced_field_length import make_balanced_field_length_problem + try: import jax except ImportError: @@ -77,12 +80,12 @@ def climb_ode(rho, S, CD0, CL0, CL_max, alpha_max, h_w, AR, e, span, T, m, v, h, return CL, q, L, D, K, F_r, v_dot, r_dot, W, v_stall, v_over_v_stall, gam_dot, h_dot -def wrap_ode_func(num_nodes, flight_mode, grad_method='jax', jax_jit=True): +def wrap_ode_func(num_nodes, mode, grad_method='jax', jax_jit=True): """ Returns the metadata from omf needed to create a new ExplciitFuncComp. """ nn = num_nodes - ode_func = runway_ode if flight_mode == 'runway' else climb_ode + ode_func = runway_ode if mode == 'runway' else climb_ode meta = (omf.wrap(ode_func) .add_input('rho', val=1.225, desc='atmospheric density at runway', units='kg/m**3') @@ -119,7 +122,7 @@ def wrap_ode_func(num_nodes, flight_mode, grad_method='jax', jax_jit=True): .add_output('v_over_v_stall', shape=(nn,), desc='stall speed ratio', units=None) ) - if flight_mode == 'runway': + if mode == 'runway': meta.add_input('mu_r', val=0.05, desc='runway friction coefficient', units=None) else: meta.add_input('gam', shape=(nn,), desc='flight path angle', units='rad') @@ -131,7 +134,10 @@ def wrap_ode_func(num_nodes, flight_mode, grad_method='jax', jax_jit=True): meta.declare_coloring('*', method=grad_method) meta.declare_partials(of='*', wrt='*', method=grad_method) - return om.ExplicitFuncComp(meta, use_jax=grad_method == 'jax', use_jit=jax_jit) + if om_version()[0] > (3, 35, 0): + return om.ExplicitFuncComp(meta, derivs_method=grad_method, use_jit=jax_jit) + else: + return om.ExplicitFuncComp(meta, use_jax=grad_method == 'jax', use_jit=jax_jit) @use_tempdirs @@ -140,181 +146,18 @@ class TestBalancedFieldFuncComp(unittest.TestCase): @unittest.skipIf(jax is None, 'requires jax and jaxlib') @require_pyoptsparse('IPOPT') def test_balanced_field_func_comp_radau(self): - self._run_problem(dm.Radau) + p = make_balanced_field_length_problem(ode_class=wrap_ode_func, + tx=dm.Radau(num_segments=3)) + + dm.run_problem(p, run_driver=True, simulate=True) + + assert_near_equal(p.get_val('traj.rto.states:r')[-1], 2197.7, tolerance=0.01) @unittest.skipIf(jax is None, 'requires jax and jaxlib') @require_pyoptsparse('IPOPT') def test_balanced_field_func_comp_gl(self): - self._run_problem(dm.GaussLobatto) - - def _run_problem(self, tx): - p = om.Problem() - - p.driver = om.pyOptSparseDriver() - p.driver.declare_coloring() - - # Use IPOPT if available, with fallback to SLSQP - p.driver.options['optimizer'] = 'IPOPT' - p.driver.options['print_results'] = True - - p.driver.opt_settings['print_level'] = 0 - p.driver.opt_settings['mu_strategy'] = 'adaptive' - - p.driver.opt_settings['bound_mult_init_method'] = 'mu-based' - p.driver.opt_settings['mu_init'] = 0.01 - p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' - - # First Phase: Brake release to V1 - both engines operable - br_to_v1 = dm.Phase(ode_class=wrap_ode_func, transcription=tx(num_segments=3), - ode_init_kwargs={'flight_mode': 'runway'}) - br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0) - br_to_v1.add_state('r', fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0) - br_to_v1.add_state('v', fix_initial=True, lower=0, ref=100.0, defect_ref=100.0) - br_to_v1.add_parameter('alpha', val=0.0, opt=False, units='deg') - br_to_v1.add_timeseries_output('*') - - # Second Phase: Rejected takeoff at V1 - no engines operable - rto = dm.Phase(ode_class=wrap_ode_func, transcription=tx(num_segments=3), - ode_init_kwargs={'flight_mode': 'runway'}) - rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) - rto.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - rto.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - rto.add_parameter('alpha', val=0.0, opt=False, units='deg') - rto.add_timeseries_output('*') - - # Third Phase: V1 to Vr - single engine operable - v1_to_vr = dm.Phase(ode_class=wrap_ode_func, transcription=tx(num_segments=3), - ode_init_kwargs={'flight_mode': 'runway'}) - v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) - v1_to_vr.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - v1_to_vr.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - v1_to_vr.add_parameter('alpha', val=0.0, opt=False, units='deg') - v1_to_vr.add_timeseries_output('*') - - # Fourth Phase: Rotate - single engine operable - rotate = dm.Phase(ode_class=wrap_ode_func, transcription=tx(num_segments=3), - ode_init_kwargs={'flight_mode': 'runway'}) - rotate.set_time_options(fix_initial=False, duration_bounds=(1.0, 5), duration_ref=1.0) - rotate.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - rotate.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - rotate.add_control('alpha', order=1, opt=True, units='deg', lower=0, upper=10, ref=10, - val=[0, 10], control_type='polynomial') - rotate.add_timeseries_output('*') - - # Fifth Phase: Climb to target speed and altitude at end of runway. - climb = dm.Phase(ode_class=wrap_ode_func, transcription=tx(num_segments=5), - ode_init_kwargs={'flight_mode': 'climb'}) - climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0) - climb.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - climb.add_state('h', fix_initial=True, lower=0, ref=1.0, defect_ref=1.0) - climb.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - climb.add_state('gam', fix_initial=True, lower=0, ref=0.05, defect_ref=0.05) - climb.add_control('alpha', opt=True, units='deg', lower=-10, upper=15, ref=10) - climb.add_timeseries_output('*') - - # Instantiate the trajectory and add phases - traj = dm.Trajectory() - p.model.add_subsystem('traj', traj) - traj.add_phase('br_to_v1', br_to_v1) - traj.add_phase('rto', rto) - traj.add_phase('v1_to_vr', v1_to_vr) - traj.add_phase('rotate', rotate) - traj.add_phase('climb', climb) - - all_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb'] - groundroll_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate'] - - # Add parameters common to multiple phases to the trajectory - traj.add_parameter('m', val=174200., opt=False, units='lbm', - desc='aircraft mass', - targets={phase: ['m'] for phase in all_phases}) - - # Handle parameters which change from phase to phase. - traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True, - desc='nominal aircraft thrust', - targets={'br_to_v1': ['T']}) - - traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_target=True, - desc='thrust under a single engine', - targets={'v1_to_vr': ['T'], 'rotate': ['T'], 'climb': ['T']}) - - traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_target=True, - desc='thrust when engines are shut down for rejected takeoff', - targets={'rto': ['T']}) - - traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_target=True, - desc='nominal runway friction coefficient', - targets={'br_to_v1': ['mu_r'], 'v1_to_vr': ['mu_r'], 'rotate': ['mu_r']}) - - traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_target=True, - desc='runway friction coefficient under braking', - targets={'rto': ['mu_r']}) - - traj.add_parameter('h_runway', val=0., opt=False, units='ft', - desc='runway altitude', - targets={phase: ['h'] for phase in groundroll_phases}) - - # Here we're omitting some constants that are common throughout all phases for the sake of brevity. - # Their correct defaults are specified in add_input calls to `wrap_ode_func`. - - # Standard "end of first phase to beginning of second phase" linkages - # Alpha changes from being a parameter in v1_to_vr to a polynomial control - # in rotate, to a dynamic control in `climb`. - traj.link_phases(['br_to_v1', 'v1_to_vr'], vars=['time', 'r', 'v']) - traj.link_phases(['v1_to_vr', 'rotate'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['rotate', 'climb'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['br_to_v1', 'rto'], vars=['time', 'r', 'v']) - - # Less common "final value of r must match at ends of two phases". - traj.add_linkage_constraint(phase_a='rto', var_a='r', loc_a='final', - phase_b='climb', var_b='r', loc_b='final', - ref=1000) - - # Define the constraints and objective for the optimal control problem - v1_to_vr.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.2, ref=100) - - rto.add_boundary_constraint('v', loc='final', equals=0., ref=100, linear=True) - - rotate.add_boundary_constraint('F_r', loc='final', equals=0, ref=100000) - - climb.add_boundary_constraint('h', loc='final', equals=35, ref=35, units='ft', linear=True) - climb.add_boundary_constraint('gam', loc='final', equals=5, ref=5, units='deg', linear=True) - climb.add_path_constraint('gam', lower=0, upper=5, ref=5, units='deg') - climb.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.25, ref=1.25) - - rto.add_objective('r', loc='final', ref=1000.0) - - # - # Setup the problem and set the initial guess - # - p.setup(check=True) - - br_to_v1.set_time_val(initial=0.0, duration=35.0) - br_to_v1.set_state_val('r', [0, 2500.0]) - br_to_v1.set_state_val('v', [0.0001, 100.0]) - br_to_v1.set_parameter_val('alpha', 0.0, units='deg') - - v1_to_vr.set_time_val(initial=35.0, duration=35.0) - v1_to_vr.set_state_val('r', [2500, 300.0]) - v1_to_vr.set_state_val('v', [100, 110.0]) - v1_to_vr.set_parameter_val('alpha', 0.0, units='deg') - - rto.set_time_val(initial=35.0, duration=1.0) - rto.set_state_val('r', [2500, 5000.0]) - rto.set_state_val('v', [110, 0.0001]) - rto.set_parameter_val('alpha', 0.0, units='deg') - - rotate.set_time_val(initial=35.0, duration=5.0) - rotate.set_state_val('r', [1750, 1800.0]) - rotate.set_state_val('v', [80, 85.0]) - rotate.set_control_val('alpha', 0.0, units='deg') - - climb.set_time_val(initial=30.0, duration=20.0) - climb.set_state_val('r', [5000, 5500.0], units='ft') - climb.set_state_val('v', [160, 170.0], units='kn') - climb.set_state_val('h', [0.0, 35.0], units='ft') - climb.set_state_val('gam', [0.0, 5.0], units='deg') - climb.set_control_val('alpha', 5.0, units='deg') + p = make_balanced_field_length_problem(ode_class=wrap_ode_func, + tx=dm.GaussLobatto(num_segments=3)) dm.run_problem(p, run_driver=True, simulate=True) diff --git a/dymos/examples/balanced_field/test/test_balanced_field_length.py b/dymos/examples/balanced_field/test/test_balanced_field_length.py index 9fc428257..16f68e1a3 100644 --- a/dymos/examples/balanced_field/test/test_balanced_field_length.py +++ b/dymos/examples/balanced_field/test/test_balanced_field_length.py @@ -10,225 +10,25 @@ import dymos as dm from dymos.examples.balanced_field.balanced_field_ode import BalancedFieldODEComp +from dymos.examples.balanced_field.balanced_field_length import make_balanced_field_length_problem from dymos.utils.misc import om_version @use_tempdirs class TestBalancedFieldLengthRestart(unittest.TestCase): - def _make_problem(self): - p = om.Problem(reports=True) - - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = 'IPOPT' - p.driver.options['invalid_desvar_behavior'] = 'ignore' - p.driver.opt_settings['print_level'] = 0 - p.driver.opt_settings['derivative_test'] = 'first-order' - - p.driver.declare_coloring() - - # First Phase: Brake release to V1 - both engines operable - br_to_v1 = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0) - br_to_v1.add_state('r', fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0) - br_to_v1.add_state('v', fix_initial=True, lower=0.0001, ref=100.0, defect_ref=100.0) - br_to_v1.add_parameter('alpha', val=0.0, opt=False, units='deg') - br_to_v1.add_timeseries_output('*') - - # Second Phase: Rejected takeoff at V1 - no engines operable - rto = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) - rto.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - rto.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0) - rto.add_parameter('alpha', val=0.0, opt=False, units='deg') - rto.add_timeseries_output('*') - - # Third Phase: V1 to Vr - single engine operable - v1_to_vr = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) - v1_to_vr.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - v1_to_vr.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0) - v1_to_vr.add_parameter('alpha', val=0.0, opt=False, units='deg') - v1_to_vr.add_timeseries_output('*') - - # Fourth Phase: Rotate - single engine operable - rotate = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - rotate.set_time_options(fix_initial=False, duration_bounds=(1.0, 5), duration_ref=1.0) - rotate.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - rotate.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0) - rotate.add_control('alpha', order=1, opt=True, units='deg', lower=0, upper=10, ref=10, - val=[0, 10], control_type='polynomial') - rotate.add_timeseries_output('*') - - # Fifth Phase: Climb to target speed and altitude at end of runway. - climb = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=5), - ode_init_kwargs={'mode': 'climb'}) - climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0) - climb.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - climb.add_state('h', fix_initial=True, lower=0.0, ref=1.0, defect_ref=1.0) - climb.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0) - climb.add_state('gam', fix_initial=True, lower=0.0, ref=0.05, defect_ref=0.05) - climb.add_control('alpha', opt=True, units='deg', lower=-10, upper=15, ref=10) - climb.add_timeseries_output('*') - - # Instantiate the trajectory and add phases - traj = dm.Trajectory() - p.model.add_subsystem('traj', traj) - traj.add_phase('br_to_v1', br_to_v1) - traj.add_phase('rto', rto) - traj.add_phase('v1_to_vr', v1_to_vr) - traj.add_phase('rotate', rotate) - traj.add_phase('climb', climb) - - # Add parameters common to multiple phases to the trajectory - traj.add_parameter('m', val=174200., opt=False, units='lbm', - desc='aircraft mass', - targets={'br_to_v1': ['m'], 'v1_to_vr': ['m'], 'rto': ['m'], - 'rotate': ['m'], 'climb': ['m']}) - - traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True, - desc='nominal aircraft thrust', - targets={'br_to_v1': ['T']}) - - traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_target=True, - desc='thrust under a single engine', - targets={'v1_to_vr': ['T'], 'rotate': ['T'], 'climb': ['T']}) - - traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_target=True, - desc='thrust when engines are shut down for rejected takeoff', - targets={'rto': ['T']}) - - traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_target=True, - desc='nominal runway friction coeffcient', - targets={'br_to_v1': ['mu_r'], 'v1_to_vr': ['mu_r'], 'rotate': ['mu_r']}) - - traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_target=True, - desc='runway friction coefficient under braking', - targets={'rto': ['mu_r']}) - - traj.add_parameter('h_runway', val=0., opt=False, units='ft', static_target=False, - desc='runway altitude', - targets={'br_to_v1': ['h'], 'v1_to_vr': ['h'], 'rto': ['h'], - 'rotate': ['h']}) - - traj.add_parameter('rho', val=1.225, opt=False, units='kg/m**3', static_target=True, - desc='atmospheric density', - targets={'br_to_v1': ['rho'], 'v1_to_vr': ['rho'], 'rto': ['rho'], - 'rotate': ['rho']}) - - traj.add_parameter('S', val=124.7, opt=False, units='m**2', static_target=True, - desc='aerodynamic reference area', - targets={'br_to_v1': ['S'], 'v1_to_vr': ['S'], 'rto': ['S'], - 'rotate': ['S'], 'climb': ['S']}) - - traj.add_parameter('CD0', val=0.03, opt=False, units=None, static_target=True, - desc='zero-lift drag coefficient', - targets={f'{phase}': ['CD0'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('AR', val=9.45, opt=False, units=None, static_target=True, - desc='wing aspect ratio', - targets={f'{phase}': ['AR'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('e', val=801, opt=False, units=None, static_target=True, - desc='Oswald span efficiency factor', - targets={f'{phase}': ['e'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('span', val=35.7, opt=False, units='m', static_target=True, - desc='wingspan', - targets={f'{phase}': ['span'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('h_w', val=1.0, opt=False, units='m', static_target=True, - desc='height of wing above CG', - targets={f'{phase}': ['h_w'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('CL0', val=0.5, opt=False, units=None, static_target=True, - desc='zero-alpha lift coefficient', - targets={f'{phase}': ['CL0'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('CL_max', val=2.0, opt=False, units=None, static_target=True, - desc='maximum lift coefficient for linear fit', - targets={f'{phase}': ['CL_max'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - # Standard "end of first phase to beginning of second phase" linkages - traj.link_phases(['br_to_v1', 'v1_to_vr'], vars=['time', 'r', 'v']) - traj.link_phases(['v1_to_vr', 'rotate'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['rotate', 'climb'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['br_to_v1', 'rto'], vars=['time', 'r', 'v']) - - # Less common "final value of r must be the match at ends of two phases". - traj.add_linkage_constraint(phase_a='rto', var_a='r', loc_a='final', - phase_b='climb', var_b='r', loc_b='final', - ref=1000) - - # Define the constraints and objective for the optimal control problem - rto.add_boundary_constraint('v', loc='final', upper=0.001, ref=100, linear=True) - - rotate.add_boundary_constraint('F_r', loc='final', equals=0, ref=100000) - - climb.add_boundary_constraint('h', loc='final', equals=35, ref=35, units='ft', linear=True) - climb.add_boundary_constraint('gam', loc='final', equals=5, ref=5, units='deg', linear=True) - climb.add_path_constraint('gam', lower=0, upper=5, ref=5, units='deg') - climb.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.2, ref=1.2) - - rto.add_objective('r', loc='final', ref=1000.0) - - # - # Setup the problem and set the initial guess - # - p.setup(check=True) - - br_to_v1.set_time_val(initial=0.0, duration=35.0) - br_to_v1.set_state_val('r', [0, 2500.0]) - br_to_v1.set_state_val('v', [0.0001, 100.0]) - br_to_v1.set_parameter_val('alpha', 0.0, units='deg') - - v1_to_vr.set_time_val(initial=35.0, duration=35.0) - v1_to_vr.set_state_val('r', [2500, 300.0]) - v1_to_vr.set_state_val('v', [100, 110.0]) - v1_to_vr.set_parameter_val('alpha', 0.0, units='deg') - - rto.set_time_val(initial=35.0, duration=1.0) - rto.set_state_val('r', [2500, 5000.0]) - rto.set_state_val('v', [110, 0.0001]) - rto.set_parameter_val('alpha', 0.0, units='deg') - - rotate.set_time_val(initial=35.0, duration=5.0) - rotate.set_state_val('r', [1750, 1800.0]) - rotate.set_state_val('v', [80, 85.0]) - rotate.set_control_val('alpha', 0.0, units='deg') - - climb.set_time_val(initial=30.0, duration=20.0) - climb.set_state_val('r', [5000, 5500.0], units='ft') - climb.set_state_val('v', [160, 170.0], units='kn') - climb.set_state_val('h', [0.0, 35.0], units='ft') - climb.set_state_val('gam', [0.0, 5.0], units='deg') - climb.set_control_val('alpha', 5.0, units='deg') - - return p - @require_pyoptsparse(optimizer='IPOPT') @unittest.skipUnless(Version(openmdao.__version__) > Version("3.23"), reason='Test requires OpenMDAO 3.23.0 or later.') def test_make_plots(self): - p = self._make_problem() + p = make_balanced_field_length_problem(ode_class=BalancedFieldODEComp, tx=dm.Radau(num_segments=3)) dm.run_problem(p, run_driver=True, simulate=True, make_plots=True) @require_pyoptsparse(optimizer='IPOPT') @unittest.skipUnless(Version(openmdao.__version__) > Version("3.23"), reason='Test requires OpenMDAO 3.23.0 or later.') def test_restart_from_sol(self): - p = self._make_problem() + p = make_balanced_field_length_problem(ode_class=BalancedFieldODEComp, tx=dm.Radau(num_segments=3)) dm.run_problem(p, run_driver=True, simulate=False) sol_db = 'dymos_solution.db' @@ -237,7 +37,7 @@ def test_restart_from_sol(self): sol_results = om.CaseReader(sol_db).get_case('final') - assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) dm.run_problem(p, run_driver=True, simulate=True, restart=sol_db) @@ -250,17 +50,17 @@ def test_restart_from_sol(self): sol_results = om.CaseReader(sol_db).get_case('final') sim_results = om.CaseReader(sim_db).get_case('final') - assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) - assert_near_equal(sim_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) + assert_near_equal(sim_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) - assert_near_equal(sol_results.get_val('traj.rto.timeseries.r')[-1], 1881, tolerance=0.01) - assert_near_equal(sim_results.get_val('traj.rto.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.rto.timeseries.r')[-1], 2197, tolerance=0.01) + assert_near_equal(sim_results.get_val('traj.rto.timeseries.r')[-1], 2197, tolerance=0.01) @require_pyoptsparse(optimizer='IPOPT') @unittest.skipUnless(Version(openmdao.__version__) > Version("3.23"), reason='Test requires OpenMDAO 3.23.0 or later.') def test_restart_from_sim(self): - p = self._make_problem() + p = make_balanced_field_length_problem(ode_class=BalancedFieldODEComp, tx=dm.Radau(num_segments=3)) dm.run_problem(p, run_driver=True, simulate=True) sol_db = 'dymos_solution.db' @@ -271,215 +71,39 @@ def test_restart_from_sim(self): sol_results = om.CaseReader(sol_db).get_case('final') - assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) dm.run_problem(p, run_driver=True, simulate=True, restart=sim_db) sol_results = om.CaseReader(sol_db).get_case('final') sim_results = om.CaseReader(sim_db).get_case('final') - assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) - assert_near_equal(sim_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) + assert_near_equal(sim_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) - assert_near_equal(sol_results.get_val('traj.rto.timeseries.r')[-1], 1881, tolerance=0.01) - assert_near_equal(sim_results.get_val('traj.rto.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.rto.timeseries.r')[-1], 2197, tolerance=0.01) + assert_near_equal(sim_results.get_val('traj.rto.timeseries.r')[-1], 2197, tolerance=0.01) @use_tempdirs class TestBalancedFieldLengthDefaultValues(unittest.TestCase): + @require_pyoptsparse(optimizer='IPOPT') def test_default_vals_stick(self): """ Make the balanced field problem without any set_val calls after setup. """ - p = om.Problem() - - # First Phase: Brake release to V1 - both engines operable - br_to_v1 = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0, - initial_val=0.0, duration_val=35.0) - br_to_v1.add_state('r', fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0, - val=br_to_v1.interp(ys=[0, 2500.0], nodes='state_input')) - br_to_v1.add_state('v', fix_initial=True, lower=0.0001, ref=100.0, defect_ref=100.0, - val=br_to_v1.interp(ys=[0.0001, 100.0], nodes='state_input')) - br_to_v1.add_parameter('alpha', val=0.0, opt=False, units='deg') - br_to_v1.add_timeseries_output('*') - - # Second Phase: Rejected takeoff at V1 - no engines operable - rto = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0, - initial_val=35.0, duration_val=35.0) - rto.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0, - val=rto.interp(ys=[2500, 5000.0], nodes='state_input')) - rto.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0, - val=rto.interp(ys=[110, 0.0001], nodes='state_input')) - rto.add_parameter('alpha', val=0.0, opt=False, units='deg') - rto.add_timeseries_output('*') - - # Third Phase: V1 to Vr - single engine operable - v1_to_vr = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0, - initial_val=35.0, duration_val=35.0) - v1_to_vr.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0, - val=v1_to_vr.interp(ys=[2500, 300.0], nodes='state_input')) - v1_to_vr.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0, - val=v1_to_vr.interp(ys=[100, 110.0], nodes='state_input')) - v1_to_vr.add_parameter('alpha', val=0.0, opt=False, units='deg') - v1_to_vr.add_timeseries_output('*') - - # Fourth Phase: Rotate - single engine operable - rotate = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - rotate.set_time_options(fix_initial=False, duration_bounds=(1.0, 5), duration_ref=1.0, - initial_val=70.0, duration_val=5.0) - rotate.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0, - val=rotate.interp(ys=[1750, 1800.0], nodes='state_input')) - rotate.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0, - val=rotate.interp(ys=[80, 85.0], nodes='state_input')) - rotate.add_control('alpha', order=1, opt=True, units='deg', lower=0, upper=10, ref=10, - val=[0, 10], control_type='polynomial') - rotate.add_timeseries_output('*') - - # Fifth Phase: Climb to target speed and altitude at end of runway. - climb = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=5), - ode_init_kwargs={'mode': 'climb'}) - climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0, - initial_val=75.0, duration_val=10.0) - climb.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0, - val=climb.interp(ys=[1800, 2500.0], nodes='state_input')) - climb.add_state('h', fix_initial=True, lower=0.0, ref=1.0, defect_ref=1.0, val=0.0) - climb.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0, - val=climb.interp(ys=[85, 90], nodes='state_input')) - climb.add_state('gam', fix_initial=True, lower=0.0, ref=0.05, defect_ref=0.05, - val=climb.interp(ys=[0, 0.05], nodes='state_input')) - climb.add_control('alpha', opt=True, units='deg', lower=-10, upper=15, ref=10, - val=climb.interp(ys=[0.01, 0.01], nodes='control_input')) - climb.add_timeseries_output('*') - - # Instantiate the trajectory and add phases - traj = dm.Trajectory() - p.model.add_subsystem('traj', traj) - traj.add_phase('br_to_v1', br_to_v1) - traj.add_phase('rto', rto) - traj.add_phase('v1_to_vr', v1_to_vr) - traj.add_phase('rotate', rotate) - traj.add_phase('climb', climb) - - # Add parameters common to multiple phases to the trajectory - traj.add_parameter('m', val=174200., opt=False, units='lbm', - desc='aircraft mass', - targets={'br_to_v1': ['m'], 'v1_to_vr': ['m'], 'rto': ['m'], - 'rotate': ['m'], 'climb': ['m']}) - - traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True, - desc='nominal aircraft thrust', targets={'br_to_v1': ['T']}) - - traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_target=True, - desc='thrust under a single engine', - targets={'v1_to_vr': ['T'], 'rotate': ['T'], 'climb': ['T']}) - - traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_target=True, - desc='thrust when engines are shut down for rejected takeoff', - targets={'rto': ['T']}) - - traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_target=True, - desc='nominal runway friction coeffcient', - targets={'br_to_v1': ['mu_r'], 'v1_to_vr': ['mu_r'], 'rotate': ['mu_r']}) - - traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_target=True, - desc='runway friction coefficient under braking', - targets={'rto': ['mu_r']}) - - traj.add_parameter('h_runway', val=0., opt=False, units='ft', static_target=False, - desc='runway altitude', - targets={'br_to_v1': ['h'], 'v1_to_vr': ['h'], 'rto': ['h'], - 'rotate': ['h']}) - - traj.add_parameter('rho', val=1.225, opt=False, units='kg/m**3', static_target=True, - desc='atmospheric density', - targets={'br_to_v1': ['rho'], 'v1_to_vr': ['rho'], 'rto': ['rho'], - 'rotate': ['rho']}) - - traj.add_parameter('S', val=124.7, opt=False, units='m**2', static_target=True, - desc='aerodynamic reference area', - targets={'br_to_v1': ['S'], 'v1_to_vr': ['S'], 'rto': ['S'], - 'rotate': ['S'], 'climb': ['S']}) - - traj.add_parameter('CD0', val=0.03, opt=False, units=None, static_target=True, - desc='zero-lift drag coefficient', - targets={f'{phase}': ['CD0'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('AR', val=9.45, opt=False, units=None, static_target=True, - desc='wing aspect ratio', - targets={f'{phase}': ['AR'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('e', val=801, opt=False, units=None, static_target=True, - desc='Oswald span efficiency factor', - targets={f'{phase}': ['e'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('span', val=35.7, opt=False, units='m', static_target=True, - desc='wingspan', - targets={f'{phase}': ['span'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('h_w', val=1.0, opt=False, units='m', static_target=True, - desc='height of wing above CG', - targets={f'{phase}': ['h_w'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('CL0', val=0.5, opt=False, units=None, static_target=True, - desc='zero-alpha lift coefficient', - targets={f'{phase}': ['CL0'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('CL_max', val=2.0, opt=False, units=None, static_target=True, - desc='maximum lift coefficient for linear fit', - targets={f'{phase}': ['CL_max'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - # Standard "end of first phase to beginning of second phase" linkages - traj.link_phases(['br_to_v1', 'v1_to_vr'], vars=['time', 'r', 'v']) - traj.link_phases(['v1_to_vr', 'rotate'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['rotate', 'climb'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['br_to_v1', 'rto'], vars=['time', 'r', 'v']) - - # Less common "final value of r must be the match at ends of two phases". - traj.add_linkage_constraint(phase_a='rto', var_a='r', loc_a='final', - phase_b='climb', var_b='r', loc_b='final', - ref=1000) - - # Define the constraints and objective for the optimal control problem - rto.add_boundary_constraint('v', loc='final', upper=0.001, ref=100, linear=True) - - rotate.add_boundary_constraint('F_r', loc='final', equals=0, ref=100000) - - climb.add_boundary_constraint('h', loc='final', equals=35, ref=35, units='ft', linear=True) - climb.add_boundary_constraint('gam', loc='final', equals=5, ref=5, units='deg', linear=True) - climb.add_path_constraint('gam', lower=0, upper=5, ref=5, units='deg') - climb.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.2, ref=1.2) - - rto.add_objective('r', loc='final', ref=1000.0) - - # - # Setup the problem and set the initial guess - # - p.setup(check=True) + p = make_balanced_field_length_problem(ode_class=BalancedFieldODEComp, tx=dm.GaussLobatto(num_segments=3)) p.run_model() - assert_near_equal(p.get_val('traj.rotate.t_initial'), 70) + assert_near_equal(p.get_val('traj.rotate.t_initial'), 35) assert_near_equal(p.get_val('traj.rotate.t_duration'), 5) - assert_near_equal(p.get_val('traj.rotate.controls:alpha'), np.array([[0, 10]]).T) - assert_near_equal(p.get_val('traj.climb.controls:alpha'), - p.model.traj.phases.climb.interp('', [0.01, 0.01], nodes='control_input')) - assert_near_equal(p.get_val('traj.climb.states:gam'), - p.model.traj.phases.climb.interp(ys=[0.0, 0.05], nodes='state_input')) - assert_near_equal(p.get_val('traj.climb.states:h'), - p.model.traj.phases.climb.interp(ys=[0.0, 0.0], nodes='state_input')) + assert_near_equal(p.get_val('traj.rotate.controls:alpha'), np.array([[0, 0]]).T) + assert_near_equal(p.get_val('traj.climb.controls:alpha', units='deg'), + p.model.traj.phases.climb.interp('', [5, 5], nodes='control_input')) + assert_near_equal(p.get_val('traj.climb.states:gam', units='deg'), + p.model.traj.phases.climb.interp(ys=[0.0, 5.0], nodes='state_input')) + assert_near_equal(p.get_val('traj.climb.states:h', units='ft'), + p.model.traj.phases.climb.interp(ys=[0.0, 35.0], nodes='state_input')) assert_near_equal(p.get_val('traj.v1_to_vr.parameters:alpha'), 0.0) diff --git a/dymos/examples/cannonball/cannonball_ode.py b/dymos/examples/cannonball/cannonball_ode.py index 7e9bded31..b9b97d1f7 100644 --- a/dymos/examples/cannonball/cannonball_ode.py +++ b/dymos/examples/cannonball/cannonball_ode.py @@ -1,14 +1,16 @@ import numpy as np -from scipy.interpolate import interp1d import openmdao.api as om +from openmdao.components.interp_util.interp import InterpND + from dymos.models.atmosphere.atmos_1976 import USatm1976Data english_to_metric_rho = om.unit_conversion('slug/ft**3', 'kg/m**3')[0] english_to_metric_h = om.unit_conversion('ft', 'm')[0] -rho_interp = interp1d(np.array(USatm1976Data.alt*english_to_metric_h, dtype=complex), - np.array(USatm1976Data.rho*english_to_metric_rho, dtype=complex), kind='linear') +rho_interp = InterpND(points=np.array(USatm1976Data.alt*english_to_metric_h), + values=np.array(USatm1976Data.rho*english_to_metric_rho), + method='slinear').interpolate class CannonballODE(om.ExplicitComponent): @@ -56,11 +58,7 @@ def compute(self, inputs, outputs): GRAVITY = 9.80665 # m/s**2 - # handle complex-step gracefully from the interpolant - if np.iscomplexobj(h): - rho = rho_interp(inputs['h']) - else: - rho = rho_interp(inputs['h']).real + rho = rho_interp(h) q = 0.5*rho*inputs['v']**2 qS = q * S diff --git a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py index 194dd79b5..ddb0f9b55 100644 --- a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py +++ b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py @@ -24,9 +24,8 @@ class TestTwoPhaseCannonballForDocs(unittest.TestCase): @require_pyoptsparse(optimizer='IPOPT') def test_two_phase_cannonball_for_docs(self): - from scipy.interpolate import interp1d - import openmdao.api as om + from openmdao.components.interp_util.interp import InterpND from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -107,16 +106,16 @@ def setup(self): self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r']) self.add_output('ke', shape=nn, units='J') - # Ask OpenMDAO to compute the partial derivatives using complex-step + # Ask OpenMDAO to compute the partial derivatives using finite-difference # with a partial coloring algorithm for improved performance, and use # a graph coloring algorithm to automatically detect the sparsity pattern. - self.declare_coloring(wrt='*', method='cs') + self.declare_coloring(wrt='*', method='fd') alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0] rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0] - self.rho_interp = interp1d(np.array(alt_data, dtype=complex), - np.array(rho_data, dtype=complex), - kind='linear') + self.rho_interp = InterpND(points=np.array(alt_data), + values=np.array(rho_data), + method='slinear') def compute(self, inputs, outputs): @@ -129,13 +128,9 @@ def compute(self, inputs, outputs): GRAVITY = 9.80665 # m/s**2 - # handle complex-step gracefully from the interpolant - if np.iscomplexobj(h): - rho = self.rho_interp(inputs['h']) - else: - rho = self.rho_interp(inputs['h']).real + rho = self.rho_interp.interpolate(h) - q = 0.5*rho*inputs['v']**2 + q = 0.5*rho*v**2 qS = q * S D = qS * CD cgam = np.cos(gam) diff --git a/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py b/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py index 96909213e..81ebd9008 100644 --- a/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py +++ b/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py @@ -9,9 +9,9 @@ class TestTwoPhaseCannonballForDocs(unittest.TestCase): @require_pyoptsparse(optimizer='SLSQP') def test_two_phase_cannonball_for_docs(self): import numpy as np - from scipy.interpolate import interp1d import openmdao.api as om + from openmdao.components.interp_util.interp import InterpND from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -92,16 +92,16 @@ def setup(self): self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r']) self.add_output('ke', shape=nn, units='J') - # Ask OpenMDAO to compute the partial derivatives using complex-step + # Ask OpenMDAO to compute the partial derivatives using finite-difference # with a partial coloring algorithm for improved performance, and use # a graph coloring algorithm to automatically detect the sparsity pattern. - self.declare_coloring(wrt='*', method='cs') + self.declare_coloring(wrt='*', method='fd') alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0] rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0] - self.rho_interp = interp1d(np.array(alt_data, dtype=complex), - np.array(rho_data, dtype=complex), - kind='linear') + self.rho_interp = InterpND(points=np.array(alt_data), + values=np.array(rho_data), + method='slinear') def compute(self, inputs, outputs): @@ -114,11 +114,7 @@ def compute(self, inputs, outputs): GRAVITY = 9.80665 # m/s**2 - # handle complex-step gracefully from the interpolant - if np.iscomplexobj(h): - rho = self.rho_interp(inputs['h']) - else: - rho = self.rho_interp(inputs['h']).real + rho = self.rho_interp.interpolate(inputs['h']) q = 0.5*rho*inputs['v']**2 qS = q * S diff --git a/dymos/examples/cannonball/test/test_load_case_missing_phase.py b/dymos/examples/cannonball/test/test_load_case_missing_phase.py index 0aa20ec74..a5aedb161 100644 --- a/dymos/examples/cannonball/test/test_load_case_missing_phase.py +++ b/dymos/examples/cannonball/test/test_load_case_missing_phase.py @@ -184,7 +184,7 @@ def test_load_case_missing_phase(self): 848.54403903, 864.60323344, 895.67083544, 926.90913594, 934.25949414, 946.77982935, 954.88116543, 955.37021704]) - assert_near_equal(p.get_val('traj.ascent.states:h').ravel(), h_expected, tolerance=1.0E-6) + assert_near_equal(p.get_val('traj.ascent.states:h').ravel(), h_expected, tolerance=1.0E-5) if __name__ == '__main__': # pragma: no cover diff --git a/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py b/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py index 1f8e3b87a..48a8d0850 100644 --- a/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py +++ b/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py @@ -3,9 +3,9 @@ from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse import numpy as np -from scipy.interpolate import interp1d import openmdao.api as om +from openmdao.components.interp_util.interp import InterpND from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -83,16 +83,16 @@ def setup(self): self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r']) self.add_output('ke', shape=nn, units='J') - # Ask OpenMDAO to compute the partial derivatives using complex-step + # Ask OpenMDAO to compute the partial derivatives using finite-difference # with a partial coloring algorithm for improved performance, and use # a graph coloring algorithm to automatically detect the sparsity pattern. - self.declare_coloring(wrt='*', method='cs') + self.declare_coloring(wrt='*', method='fd') alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0] rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0] - self.rho_interp = interp1d(np.array(alt_data, dtype=complex), - np.array(rho_data, dtype=complex), - kind='linear', bounds_error=False, fill_value='extrapolate') + self.rho_interp = InterpND(points=np.array(alt_data), + values=np.array(rho_data), + method='slinear', extrapolate=True) def compute(self, inputs, outputs): @@ -105,11 +105,7 @@ def compute(self, inputs, outputs): GRAVITY = 9.80665 # m/s**2 - # handle complex-step gracefully from the interpolant - if np.iscomplexobj(h): - rho = self.rho_interp(inputs['h']) - else: - rho = self.rho_interp(inputs['h']).real + rho = self.rho_interp.interpolate(inputs['h']) q = 0.5 * rho * inputs['v'] ** 2 qS = q * S diff --git a/pyproject.toml b/pyproject.toml index a7c048d95..d75b1c4c9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,6 +33,7 @@ all = [ "dymos[docs]", "dymos[notebooks]", "dymos[test]", + "dymos[numba]", ] docs = [ "bokeh", @@ -62,6 +63,8 @@ test = [ "playwright>=1.20", "pycodestyle", "testflo>=1.3.6", +] +numba = [ "numba", ]