[petsc-users] Asking examples about solving DAE in python

[Xiong, Jing]

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 <jxiong at anl.gov>
Sent: Monday, January 24, 2022 10:03 PM
To: petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
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 <hongzhang at anl.gov>
Sent: Thursday, January 20, 2022 3:05 PM
To: Xiong, Jing <jxiong at anl.gov>
Cc: petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>; Zhao, Dongbo <dongbo.zhao at anl.gov>; Hong, Tianqi <thong at anl.gov>
Subject: Re: [petsc-users] Asking examples about solving DAE in python



On Jan 20, 2022, at 4:13 PM, Xiong, Jing via petsc-users <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>> 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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220125/00a7a168/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 30250 bytes
Desc: image.png
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220125/00a7a168/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 21718 bytes
Desc: image.png
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220125/00a7a168/attachment-0003.png>