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

[Zhang, Hong]

Oops. Some typos in the index. They are corrected in below.

        x[0] = x0
        x[1] = y0
        x[2] = 0
        x[3] = 0
        x[4] = lambdaval

Note that these initial values are randomly picked. Ideally, you should make sure the initial condition is consistent for DAEs. In other words, the algebraic equations (the last equation in your case) should be satisfied when you choose initial values.

Hong

> On Jan 25, 2022, at 2:29 PM, Zhang, Hong via petsc-users <petsc-users at mcs.anl.gov> wrote:
> 
> 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> 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:
>> 	• How to set the initial value for xdot?
>> 	• 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>
>> 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> 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:
>>> 	• Is petsc4py the right package to use?
>> 
>> Yes, you need petsc4py for python. 
>> 
>>> 	• 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/.
>> 	• How to solve an ODE with explicit methods.
>> 	• How to solve an ODE/DAE with implicit methods.
>> 	• How to use TSAdjoint to calculate adjoint sensitivities.
>> 	• 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). 
>> 
>>> 	• 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
>