Sorry for the confusion. I was implemented following the pendulum example given by https://github.com/bmcage/odes/blob/master/ipython_examples/Planar%20Pendulum%20as%20DAE.ipynb.
The formula is [cid:b49e0bc9-f98e-44a8-a55a-caa0fd0522e3] ________________________________ From: Xiong, Jing <[email protected]> Sent: Monday, January 24, 2022 10:03 PM To: [email protected] <[email protected]> Subject: Re: [petsc-users] Asking examples about solving DAE in python Good morning, Thanks for all your help. I'm trying to implement a toy example solving a DAE for a Pendulum system and I got two questions: 1. How to set the initial value for xdot? 2. I got the following error information: * Assertion failed: (!PyErr_Occurred()), function __Pyx_PyCFunction_FastCall, file src/petsc4py.PETSc.c, line 359099. The system equation is given in https://en.wikipedia.org/wiki/Differential-algebraic_system_of_equations The formula I used is shown as follows: [cid:3001ee1d-031e-4500-b5e0-bc14e99812b5] I also attached my code below: import sys, petsc4py petsc4py.init(sys.argv) import numpy as np from petsc4py import PETSc class Pendulum(object): n = 5 comm = PETSc.COMM_SELF def initialCondition(self, x): # mu = self.mu_ l = 1.0 m = 1.0 g = 1.0 #initial condition theta0= np.pi/3 #starting angle x0=np.sin(theta0) y0=-(l-x0**2)**.5 lambdaval = 0.1 x[0] = x0 x[1] = y0 x[4] = lambdaval x.assemble() def evalIFunction(self, ts, t, x, xdot, f): f.setArray ([x[2]-xdot[0]], [x[3]-xdot[1]], [-xdot[2]+x[4]*x[0]/m], [-xdot[3]+x[4]*x[1]/m-g], [x[2]**2 + x[3]**2 + (x[0]**2 + x[1]**2)/m*x[4] - x[1] * g]) OptDB = PETSc.Options() ode = Pendulum() x = PETSc.Vec().createSeq(ode.n, comm=ode.comm) f = x.duplicate() ts = PETSc.TS().create(comm=ode.comm) ts.setType(ts.Type.CN) ts.setIFunction(ode.evalIFunction, f) ts.setSaveTrajectory() ts.setTime(0.0) ts.setTimeStep(0.001) ts.setMaxTime(0.5) ts.setMaxSteps(1000) ts.setExactFinalTime(PETSc.TS.ExactFinalTime.MATCHSTEP) ts.setFromOptions() ode.initialCondition(x) ts.solve(x) Best, Jing ________________________________ From: Zhang, Hong <[email protected]> Sent: Thursday, January 20, 2022 3:05 PM To: Xiong, Jing <[email protected]> Cc: [email protected] <[email protected]>; Zhao, Dongbo <[email protected]>; Hong, Tianqi <[email protected]> Subject: Re: [petsc-users] Asking examples about solving DAE in python On Jan 20, 2022, at 4:13 PM, Xiong, Jing via petsc-users <[email protected]<mailto:[email protected]>> wrote: Hi, I hope you are well. 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. The only example I got right now is for solving ODEs in the paper: PETSc/TS: A Modern Scalable ODE/DAE Solver Library. And I got the following questions: 1. Is petsc4py the right package to use? Yes, you need petsc4py for python. 1. Could you give an example for solving DAEs in python? 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/. 1. How to solve an ODE with explicit methods. 2. How to solve an ODE/DAE with implicit methods. 3. How to use TSAdjoint to calculate adjoint sensitivities. 4. 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). 1. Is Jacobian must be specified? If not, could your show an example for solving DAEs without specified Jacobian in python? 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. Hong Thank you for your help. Best, Jing
