[petsc-users] Question about how to solve the DAE step by step using PETSC4PY.

Tue Apr 5 22:15:43 CDT 2022

Hi Jing,

Please send your message to the mailing list when reply, so our response time will be fastest and the messages could potentially be helpful to other users.

1. You need to set up the final time step properly with ts.setExactFinalTime(PETSc.TS.ExactFinalTime.MATCHSTEP) or the command line option -ts_exact_final_time match step (need to add ts.setFromOptions() to enable options at run time). More details can be found in https://petsc.org/release/docs/manualpages/TS/TSSetExactFinalTime.html . In method 2, TS is actually using the interpolate option, which means TS may step over the final time you specify and use interpolation to generate the solution at the specified time.

2. It seems that you want to do a single step from t1 to t2 in each iteration. This can be achieved in several different ways. One way is to add ts.setTimeStep(0.0002) in the loop. It is fine to use any step size larger than 0.0001 because TS would adapt the step size to match the final time t2. But do not set it to 0.0001 since the round-off error may result in an additional tiny step in order to match t2 exactly. Another way is to add ts.setTimeStep(index) before ts.solve() in the loop. This controls the max number of steps incrementally so that only one time step is taken in each iteration.

After these fixes, method 2 will give the same results as method 1.

Hong (Mr.)

> On Apr 5, 2022, at 8:27 PM, Xiong, Jing <jxiong at anl.gov> wrote:
> 
> Hello Dr. Zhang,
>  
> Thank you for your reply.
> I have implemented an toy example using traditional way (method 1 in code: not step by step) and step by step solve (method 2 in code). However, I couldn’t get same results. The code and the simulation results are attached below.  Could you help me with this?
>  
> Best,
> Jing
> 
> from __future__ import print_function, division
> import sys, petsc4py
> petsc4py.init(sys.argv)
> import numpy as np
> import matplotlib.pyplot as plt
> from petsc4py import PETSc
> 
> #
> class twobus3phase(object):
>     def __init__(self,v):
>         self.n = 18
>         self.comm = PETSc.COMM_SELF
>         self.v = v
> 
>     def evalIFunction(self, ts, t, x, xdot, result):
>         result[0] = self.v * np.cos(2*np.pi*60*t) - x[0] - (R[0][0] * x[3] + R[0][1] * x[4] + R[0][2] * x[5]) - (L[0][0] * xdot[3] + L[0][1] * xdot[4] + L[0][2] * xdot[5])
>         result[1] = self.v * np.cos(2*np.pi*60*t - 2*np.pi/3) - x[1] - (R[1][0] * x[3] + R[1][1] * x[4] + R[1][2] * x[5]) - (L[1][0] * xdot[3] + L[1][1] * xdot[4] + L[1][2] * xdot[5])
>         result[2] = self.v * np.cos(2*np.pi*60*t + 2*np.pi/3) - x[2] - (R[2][0] * x[3] + R[2][1] * x[4] + R[2][2] * x[5]) - (L[2][0] * xdot[3] + L[2][1] * xdot[4] + L[2][2] * xdot[5])
>         result[3] = x[12] - (C[0][0] * xdot[0] + C[0][1] * xdot[1] + C[0][2] * xdot[2])
>         result[4] = x[13] - (C[1][0] * xdot[0] + C[1][1] * xdot[1] + C[1][2] * xdot[2])
>         result[5] = x[14] - (C[2][0] * xdot[0] + C[2][1] * xdot[1] + C[2][2] * xdot[2])
>         result[6] = x[0] - loadVec[1][1] * xdot[6]
>         result[7] = x[1] - loadVec[1][1] * xdot[7]
>         result[8] = x[2] - loadVec[1][1] * xdot[8]
>         result[9] = x[0] - loadVec[1][0] * x[9]
>         result[10] = x[1] - loadVec[1][0] * x[10]
>         result[11] = x[2] - loadVec[1][0] * x[11]
>         result[12] = x[12] - x[3] + x[15] + inputVec[1][0]
>         result[13] = x[13] - x[4] + x[16] + inputVec[1][1]
>         result[14] = x[14] - x[5] + x[17] + inputVec[1][2]
>         result[15] = x[15] - x[6] - x[9]
>         result[16] = x[16] - x[7] - x[10]
>         result[17] = x[17] - x[8] - x[11]
>         result.assemble()
> 
> def RLCEqua(R0, R1):
>     return [[(2 * R1 + R0) / 3, (R0 - R1) / 3, (R0 - R1) / 3],
>             [(R0 - R1) / 3, (2 * R1 + R0) / 3, (R0 - R1) / 3],
>             [(R0 - R1) / 3, (R0 - R1) / 3, (2 * R1 + R0) / 3]
>             ]
> 
> def getC(i):
>     C = [[0, 0, 0]] * 3
>     for j in range(0, len(lineMat[0])):
>         if lineMat[i][j]:
>             C0, C1 = CValue[lineMat[i][j]-1]
>             C += RLCEqua(C0, C1)
>     return C
> #
> Sbase = 10
> Vbase = 230 / np.sqrt(3) * np.sqrt(2)
> Zbase = Vbase * Vbase / Sbase
> 
> lineMat = [[0,1],[1,0]]
> loadVec = [[],[5.2454e3, 27.83]]
> inputVec = [[],[0,0,0]]
> RValue = [[47.22, 31.02]]
> LValue = [[0.0799, 0.2444 ]]
> CValue = [[2.961e-7, 8.883e-7 ]]
> 
> num = 1
> for i in range(len(lineMat)):
>     for j in range(i, len(lineMat[0])):
>         if lineMat[i][j]:
>             lineMat[i][j] = num
>             lineMat[j][i] = num
>             num += 1
> 
> def perUnit(list, base):
>     for x in list:
>         if x:
>             for i in range(len(x)):
>                 x[i] = x[i]/base
>     return list
> 
> loadVec = perUnit(loadVec, Zbase)
> RValue = perUnit(RValue, Zbase)
> LValue = perUnit(LValue, Zbase)
> CValue = perUnit(CValue, 1/Zbase)
> 
> R0, R1 = RValue[0]
> L0, L1 = LValue[0]
> R = RLCEqua(R0, R1)
> L = RLCEqua(L0, L1)
> C = getC(1)
> 
> if Vbase == 1:
>     v = 1.06 * 230 / np.sqrt(3) * np.sqrt(2)
> else:
>     v = 1.06
> 
> OptDB = PETSc.Options()
> 
> ode = twobus3phase(v)
> 
> x = PETSc.Vec().createSeq(ode.n, comm=ode.comm)
> f = x.duplicate()
> 
> ts = PETSc.TS().create(comm=ode.comm)
> ts.setIFunction(ode.evalIFunction, f)
> 
> # ## method 1:
> # history = []
> # def monitor(ts , i, t, x):
> #     xx = x[:]. tolist ()
> #     history.append ((i, t, xx))
> # ts.setMonitor(monitor)
> # ts.setTime (0.0)
> # ts.setTimeStep (0.0001)
> # ts.setMaxTime (0.1)
> #
> 
> ## method 2:
> stop_t  = np.linspace(0,0.1,1000)
> lastind = len(stop_t)
> timeInit = 0.
> history = []
> for index,time in enumerate(stop_t):
>     if index == 0:
>         continue
> 
>     t1 = timeInit
>     t2 = time
>     ts.setTime(t1)
>     ts.setMaxTime(t2)
>     ts.solve(x)
> 
>     xx = x[:].tolist()
>     history.append((index, time, xx))
> 
>     timeInit = time
> 
> ## plot
> tt = np.asarray ([v[1] for v in history ])
> xx = np.asarray ([v[2] for v in history ])
> 
> plt.figure()
> plt.plot(tt, xx[:,0], label='v')
> plt.plot(tt, xx[:,6], label='load_il')
> plt.plot(tt, xx[:,9], label='load_ir')
> plt.legend()
> plt.show()
>  
> Method 1 result:
> <image001.png>
>  
> Method 2 result:
>  
> <image002.png>
>  
> 
> From: Zhang, Hong <hongzhang at anl.gov>
> Date: Sunday, March 27, 2022 at 7:19 PM
> To: Xiong, Jing <jxiong at anl.gov>
> Cc: petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
> Subject: Re: [petsc-users] Question about how to solve the DAE step by step using PETSC4PY.
> 
>  
> 
> 
> On Mar 25, 2022, at 1:14 PM, Xiong, Jing via petsc-users <petsc-users at mcs.anl.gov> wrote:
>  
> Good afternoon,
>  
> Thanks for all your help.
> I got a question about how to solve the DAE step by step using PETSC4PY.
> I got two DAE systems, let's say f1(x1, u1) and f2(x2, u2). The simplified algorithm I need to implement is as follows:
> Step 1: solve f1 for 1 step.
> Step 2: use part of x1 as input for f2, which means u2 = x1[:n].
> Step 3: solve f2 for one step with u2 = x1[:n].
> Step 4: use part of x2 as input for f1, which means u1 = x2[:m].
>  
> I'm not able to find any examples of how to use PETSC4PY in such scenarios. If using the "scikits.odes.dae" package, it is like:
> daesolver.init_step(timeInit, XInit, XpInit)
> daesolver.step(time)
>  
> Jing,
>  
> You can certainly do the same thing in petsc4py with
>  
> ts.setMaxTime(t1)
> ts.solve() // integrate over time interval [0,t1]
> ts.setTime(t1)
> ts.setMaxTime(t2)
> ts.solve()  // integrate over time interval [t1,t2]
> ...
>  
> There are also APIs that allow you to control the step size (ts.setTimeStep) and the final step behavior (ts.setExactFinalTime).
> The C example src/ts/tutorials/power_grid/stability_9bus/ex9busdmnetwork.c might be helpful to you. Most of the C examples can be reproduced in python with similar APIs.
>  
> Hong
> 
> 
> Please let me know if there are any examples I can refer to. Thank you.
>  
> Best,
> Jing