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
> <[email protected]> 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
>> <[email protected]> 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 <[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]> 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
>
