<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<style type="text/css" style="display:none;"> P {margin-top:0;margin-bottom:0;} </style>
</head>
<body dir="ltr">
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
Sorry for the confusion. I was implemented following the pendulum example given by <a href="https://github.com/bmcage/odes/blob/master/ipython_examples/Planar%20Pendulum%20as%20DAE.ipynb" id="LPNoLPOWALinkPreview">https://github.com/bmcage/odes/blob/master/ipython_examples/Planar%20Pendulum%20as%20DAE.ipynb</a>.</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
</div>
<div class="_Entity _EType_OWALinkPreview _EId_OWALinkPreview _EReadonly_1"></div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
The formula is</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
<img style="width: 211px; height: 158px; max-width: initial;" class="w-384 h-288" width="211" height="158" size="21718" contenttype="image/png" data-outlook-trace="F:1|T:1" src="cid:b49e0bc9-f98e-44a8-a55a-caa0fd0522e3"><br>
</div>
<div id="appendonsend"></div>
<hr style="display:inline-block;width:98%" tabindex="-1">
<div id="divRplyFwdMsg" dir="ltr"><font face="Calibri, sans-serif" style="font-size:11pt" color="#000000"><b>From:</b> Xiong, Jing <jxiong@anl.gov><br>
<b>Sent:</b> Monday, January 24, 2022 10:03 PM<br>
<b>To:</b> petsc-users@mcs.anl.gov <petsc-users@mcs.anl.gov><br>
<b>Subject:</b> Re: [petsc-users] Asking examples about solving DAE in python</font>
<div> </div>
</div>
<style type="text/css" style="display:none">
<!--
p
{margin-top:0;
margin-bottom:0}
-->
</style>
<div dir="ltr">
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
Good morning,</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
<br>
</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
Thanks for all your help.</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
I'm trying to implement a toy example solving a DAE for a Pendulum system and I got two questions:</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
<ol>
<li><span>How to set the initial value for xdot?</span></li><li><span style="color:rgb(0,0,0); font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt">I got the following error information:</span><br>
</li></ol>
</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
<ul>
<li><span><b style="font-size:12pt; font-style:inherit; font-variant-ligatures:inherit; font-variant-caps:inherit">Assertion failed: (!PyErr_Occurred()), function __Pyx_PyCFunction_FastCall, file src/petsc4py.PETSc.c, line 359099.</b><br>
</span></li></ul>
</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
<span style="color:rgb(0,0,0); font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt">The system equation is given in </span><a href="https://en.wikipedia.org/wiki/Differential-algebraic_system_of_equations" id="LPNoLPOWALinkPreview" style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt">https://en.wikipedia.org/wiki/Differential-algebraic_system_of_equations</a><br>
</div>
<div class="x__Entity x__EType_OWALinkPreview x__EId_OWALinkPreview x__EReadonly_1">
</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
The formula I used is shown as follows:</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
<img class="x_w-542 x_h-320" width="243" height="144" size="30250" style="width:243.9px; height:144px; max-width:initial" data-outlook-trace="F:2|T:2" src="cid:3001ee1d-031e-4500-b5e0-bc14e99812b5"><br>
</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
<br>
</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
I also attached my code below:</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
<pre style="background-color:#ffffff; color:#080808; font-family:'JetBrains Mono',monospace; font-size:9.8pt"><span style="color:#0033b3">import </span>sys, petsc4py<br>petsc4py.init(sys.argv)<br><span style="color:#0033b3">import </span>numpy <span style="color:#0033b3">as </span>np<br><br><span style="color:#0033b3">from </span>petsc4py <span style="color:#0033b3">import </span>PETSc<br><br><span style="color:#0033b3">class </span><span style="color:#000000">Pendulum</span>(object):<br> n = <span style="color:#1750eb">5<br></span><span style="color:#1750eb"> </span>comm = PETSc.COMM_SELF<br> <span style="color:#0033b3">def </span><span style="color:#00627a">initialCondition</span>(<span style="color:#94558d">self</span>, x):<br> <span style="color:#8c8c8c; font-style:italic"># mu = self.mu_<br></span><span style="color:#8c8c8c; font-style:italic"> </span>l = <span style="color:#1750eb">1.0<br></span><span style="color:#1750eb"> </span><span style="color:#808080">m </span>= <span style="color:#1750eb">1.0<br></span><span style="color:#1750eb"> </span><span style="color:#808080">g </span>= <span style="color:#1750eb">1.0<br></span><span style="color:#1750eb"> </span><span style="color:#8c8c8c; font-style:italic">#initial condition<br></span><span style="color:#8c8c8c; font-style:italic"> </span>theta0= np.pi/<span style="color:#1750eb">3 </span><span style="color:#8c8c8c; font-style:italic">#starting angle<br></span><span style="color:#8c8c8c; font-style:italic"> </span>x0=np.sin(theta0)<br> y0=-(l-x0**<span style="color:#1750eb">2</span>)**<span style="color:#1750eb">.5<br></span><span style="color:#1750eb"> </span>lambdaval = <span style="color:#1750eb">0.1<br></span><span style="color:#1750eb"> </span>x[<span style="color:#1750eb">0</span>] = x0<br> x[<span style="color:#1750eb">1</span>] = y0<br> x[<span style="color:#1750eb">4</span>] = lambdaval<br> x.assemble()<br><br> <span style="color:#0033b3">def </span><span style="color:#00627a">evalIFunction</span>(<span style="color:#94558d">self</span>, <span style="color:#808080">ts</span>, <span style="color:#808080">t</span>, x, xdot, f):<br> f.setArray ([x[<span style="color:#1750eb">2</span>]-xdot[<span style="color:#1750eb">0</span>]],<br> [x[<span style="color:#1750eb">3</span>]-xdot[<span style="color:#1750eb">1</span>]],<br> [-xdot[<span style="color:#1750eb">2</span>]+x[<span style="color:#1750eb">4</span>]*x[<span style="color:#1750eb">0</span>]/m],<br> [-xdot[<span style="color:#1750eb">3</span>]+x[<span style="color:#1750eb">4</span>]*x[<span style="color:#1750eb">1</span>]/m-g],<br> [x[<span style="color:#1750eb">2</span>]**<span style="color:#1750eb">2 </span>+ x[<span style="color:#1750eb">3</span>]**<span style="color:#1750eb">2 </span>+ (x[<span style="color:#1750eb">0</span>]**<span style="color:#1750eb">2 </span>+ x[<span style="color:#1750eb">1</span>]**<span style="color:#1750eb">2</span>)/m*x[<span style="color:#1750eb">4</span>] - x[<span style="color:#1750eb">1</span>] * g])<br><br>OptDB = PETSc.Options()<br>ode = Pendulum()<br><br>x = PETSc.Vec().createSeq(ode.n, <span style="color:#660099">comm</span>=ode.comm)<br>f = x.duplicate()<br><br>ts = PETSc.TS().create(<span style="color:#660099">comm</span>=ode.comm)<span style="color:#8c8c8c; font-style:italic"><br></span><span style="color:#8c8c8c; font-style:italic"><br></span>ts.setType(ts.Type.CN)<br>ts.setIFunction(ode.evalIFunction, f)<br><br>ts.setSaveTrajectory()<br>ts.setTime(<span style="color:#1750eb">0.0</span>)<br>ts.setTimeStep(<span style="color:#1750eb">0.001</span>)<br>ts.setMaxTime(<span style="color:#1750eb">0.5</span>)<br>ts.setMaxSteps(<span style="color:#1750eb">1000</span>)<br>ts.setExactFinalTime(PETSc.TS.ExactFinalTime.MATCHSTEP)<br><span style="color:#8c8c8c; font-style:italic"><br></span>ts.setFromOptions()<br>ode.initialCondition(x)<br>ts.solve(x)</pre>
<br>
</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
Best,<br>
Jing</div>
<div id="x_appendonsend"></div>
<hr tabindex="-1" style="display:inline-block; width:98%">
<div id="x_divRplyFwdMsg" dir="ltr"><font face="Calibri, sans-serif" color="#000000" style="font-size:11pt"><b>From:</b> Zhang, Hong <hongzhang@anl.gov><br>
<b>Sent:</b> Thursday, January 20, 2022 3:05 PM<br>
<b>To:</b> Xiong, Jing <jxiong@anl.gov><br>
<b>Cc:</b> petsc-users@mcs.anl.gov <petsc-users@mcs.anl.gov>; Zhao, Dongbo <dongbo.zhao@anl.gov>; Hong, Tianqi <thong@anl.gov><br>
<b>Subject:</b> Re: [petsc-users] Asking examples about solving DAE in python</font>
<div> </div>
</div>
<div class="" style="word-wrap:break-word; line-break:after-white-space"><br class="">
<div><br class="">
<blockquote type="cite" class="">
<div class="">On Jan 20, 2022, at 4:13 PM, Xiong, Jing via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" class="">petsc-users@mcs.anl.gov</a>> wrote:</div>
<br class="x_x_Apple-interchange-newline">
<div class="">
<div class="" style="font-style:normal; font-variant-caps:normal; font-weight:normal; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none; font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt">
<div class="" style="margin-top:0pt; margin-bottom:0pt; color:rgb(14,16,26); background-color:transparent">
<span class="" style="background-color:transparent">Hi,</span></div>
<div class="" style="margin-top:0pt; margin-bottom:0pt; color:rgb(14,16,26); background-color:transparent">
<br class="">
</div>
<div class="" style="margin-top:0pt; margin-bottom:0pt; color:rgb(14,16,26); background-color:transparent">
<span class="" style="background-color:transparent">I hope you are well.</span></div>
<div class="" style="margin-top:0pt; margin-bottom:0pt; color:rgb(14,16,26); background-color:transparent">
<span class="" style="background-color:transparent">I'm very interested in PETSc and want to explore the possibility of whether it could solve Differential-algebraic equations (DAE) in python. I know there are great examples in C, but I'm struggling to connect
the examples in python. </span></div>
<div class="" style="margin-top:0pt; margin-bottom:0pt; color:rgb(14,16,26); background-color:transparent">
<br class="">
</div>
<div class="" style="margin-top:0pt; margin-bottom:0pt; color:rgb(14,16,26); background-color:transparent">
<span class="" style="background-color:transparent">The only example I got right now is for solving ODEs in the paper: PETSc/TS: A Modern Scalable ODE/DAE Solver Library. </span></div>
<div class="" style="margin-top:0pt; margin-bottom:0pt; color:rgb(14,16,26); background-color:transparent">
<span class="" style="background-color:transparent">And I got the following questions:</span></div>
<ol class="" style="color:rgb(14,16,26); background-color:transparent; margin-top:0pt; margin-bottom:0pt">
<li class="" style="background-color:transparent; margin-top:0pt; margin-bottom:0pt; list-style-type:decimal">
<span class="" style="background-color:transparent; margin-top:0pt; margin-bottom:0pt"><b class="">Is petsc4py the right package to use?</b></span></li></ol>
</div>
</div>
</blockquote>
<div><br class="">
</div>
<div>Yes, you need petsc4py for python. </div>
<br class="">
<blockquote type="cite" class="">
<div class="">
<div class="" style="font-style:normal; font-variant-caps:normal; font-weight:normal; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none; font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt">
<ol class="" start="2" style="color:rgb(14,16,26); background-color:transparent; margin-top:0pt; margin-bottom:0pt">
<li class="" style="background-color:transparent; margin-top:0pt; margin-bottom:0pt; list-style-type:decimal">
<span class="" style="background-color:transparent; margin-top:0pt; margin-bottom:0pt"><b class="">Could you give an example for solving DAEs in python?</b></span></li></ol>
</div>
</div>
</blockquote>
<div><br class="">
</div>
<div>src/binding/petsc4py/demo/ode/vanderpol.py gives a simple example to demonstrate a variety of PETSc capabilities. The corresponding C version of this example is ex20adj.c in src/ts/tutorials/.</div>
<div>
<ol class="x_x_MailOutline">
<li class="">How to solve an ODE with explicit methods.</li><li class="">How to solve an ODE/DAE with implicit methods.</li><li class="">How to use TSAdjoint to calculate adjoint sensitivities.</li><li class="">How to do a manual matrix-free implementation (e.g. you may already have a way to differentiate your RHS function to generate the Jacobian-vector product). </li></ol>
</div>
<br class="">
<blockquote type="cite" class="">
<div class="">
<div class="" style="font-style:normal; font-variant-caps:normal; font-weight:normal; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none; font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt">
<ol class="" start="3" style="color:rgb(14,16,26); background-color:transparent; margin-top:0pt; margin-bottom:0pt">
<li class="" style="background-color:transparent; margin-top:0pt; margin-bottom:0pt; list-style-type:decimal">
<span class="" style="background-color:transparent; margin-top:0pt; margin-bottom:0pt"><b class="">Is Jacobian must be specified? If not, could your show an example for solving DAEs without specified Jacobian in python?</b></span></li></ol>
</div>
</div>
</blockquote>
<div><br class="">
</div>
PETSc can do finite-difference approximations to generate the Jacobian matrix automatically. This may work nicely for small-sized problems. You can also try to use an AD tool to produce the Jacobian-vector product and use it in a matrix-free implementation.</div>
<div><br class="">
</div>
<div>Hong</div>
<div><br class="">
<blockquote type="cite" class="">
<div class="">
<div class="" style="font-style:normal; font-variant-caps:normal; font-weight:normal; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none; font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt">
<div class=""><span class="" style="background-color:transparent; margin-top:0pt; margin-bottom:0pt"><br class="">
</span></div>
<div class=""><span class="" style="background-color:transparent; margin-top:0pt; margin-bottom:0pt">Thank you for your help.</span></div>
<div class=""><span class="" style="background-color:transparent; margin-top:0pt; margin-bottom:0pt"><b class=""><br class="">
</b></span></div>
<div class=""><span class="" style="background-color:transparent; margin-top:0pt; margin-bottom:0pt">Best,</span></div>
<div class=""><span class="" style="background-color:transparent; margin-top:0pt; margin-bottom:0pt">Jing</span></div>
</div>
</div>
</blockquote>
</div>
<br class="">
</div>
</div>
</body>
</html>