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

[Zhang, Hong]

The following code should work.

class Pendulum(object):
    n = 5
    comm = PETSc.COMM_SELF
    l = 1.0
    m = 1.0
    g = 1.0
    def initialCondition(self, x):
        # mu = self.mu_
        #initial condition
        theta0= np.pi/3 #starting angle
        x0=np.sin(theta0)
        y0=-(self.l-x0**2)**.5
        lambdaval = 0.1
        x[0] = x0
        x[1] = y0
        x[2] = 0
        x[4] = 0
        x[5] = lambdaval
        x.assemble()

    def evalIFunction(self, ts, t, x, xdot, f):
        f[0] = x[2]-xdot[0]
        f[1] = x[3]-xdot[1]
        f[2] = -xdot[2]+x[4]*x[0]/self.m
        f[3] = -xdot[3]+x[4]*x[1]/self.m-self.g
        f[4] = [x[2]**2 + x[3]**2 + (x[0]**2 + x[1]**2)/self.m*x[4] - x[1] * self.g]
        f.assemble()

You do not need to initialize xdot. Only x needs to be initialized.

Hong

On Jan 25, 2022, at 12:03 AM, Xiong, Jing via petsc-users <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>> wrote:

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:
<image.png>

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<mailto:hongzhang at anl.gov>>
Sent: Thursday, January 20, 2022 3:05 PM
To: Xiong, Jing <jxiong at anl.gov<mailto:jxiong at anl.gov>>
Cc: petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov> <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>>; Zhao, Dongbo <dongbo.zhao at anl.gov<mailto:dongbo.zhao at anl.gov>>; Hong, Tianqi <thong at anl.gov<mailto: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/02c746ad/attachment-0001.html>