Iteration time

[retter-berkeley]

I am using py-pde in a non-standard way to try and use it to explore (through RL) controllers for PDE/distributed parameter systems. To allow the controller to act on the system, I am looping over the pde with a short time horizon, then applying the control inputs to the resulting state, and using that state as the basis for another short pde run, as in the following pseudo code:

while error > epsilon:
   use current state to determine control input
   state = state + control input
   solve pd (eq.solve) on updated state, time range = 0.2, adaptive time step, no tracker
   update state with new pde solution
   compare current state with desired state

The problem is the amount of time it takes py-pde to run eq.solve (currently using Diffusion) grows (like 5-10x) as I go through thousands of these iterations. I cannot figure out why. I start with taking about 30 sec to go through 50 iterations and after a couple hundred of these I am up to around 100-120 sec to complete 50 iterations. I need such high iterations for training reinforcement learning. I have tried garbage collect, and deleting variables like grid, state, etc. I know the time . The only way I have found to get back down to 30s iteration time is to restart the kernel in my notebook.

I know py-pde wasn't designed to be used (abused) like this, but is there something I can do to keep iteration time down? Is there some variable I can delete, or a memory leak somewhere? I can't figure out why the iteration time grows, or a more workflow-friendly way to restore it.

[david-zwicker]

Thanks for reporting – this sounds like a memory leak due to caching, potential inside numba. You might be able to confirm this if you watch the memory consumption of the process.

In any case, there might be a better work-around, which likely will also lead to better performance: You could create a stepper function directly:

grid = UnitGrid([32, 32])  # generate grid
state = ScalarField.random_uniform(grid, 0.2, 0.3)  # generate initial condition
eq = pde.DiffusionPDE()
solver = pde.ExplicitSolver(eq, scheme="euler", adaptive=True)  # initialize a solver
stepper = solver.make_stepper(state, dt=1e-3)  # initial time step

The stepper can then be used inside your loop:

while error > epsilon:
   use current state to determine control input
   state = state + control input
   t_final = stepper(state, t_current, t_current+0.2)
   compare current state with desired state

Note that the stepper function only needs to be created once. The stepper updates the state directly, so if you need the state before the stepper, you need to copy it. The stepper returns the final time, which might deviate slightly from the requested end time, particularly if fixed time steps are used.

[asmahr]

My question fits well to this post because I am running into a memory issue when I am solving the pdes in a loop. You can observe this behavior already for this simplified code

grid = UnitGrid([32, 32])
state = ScalarField.random_uniform(grid, 0.2, 0.3)  
eq = DiffusionPDE() 

for i in range(0,5):
     eq.solve(state,t_range=(1e2, 1e5),dt=0.1,adaptive=False, tracker=None, solver="explicit",ret_info = False)
    gc.collect()

As you mentioned before, I also assume some cache isn’t being freed. How can I solve the issue? For me, the main benefit of using eq.solve() is that it has the tracker argument. For my calculation, this argument is essential because it saves the states at specific times during the calculation to an external file.

In addition, reading the documentation it is not completely clear why the stepper function can speed up calculations.

[asmahr]

Sure, in my code, I use a tracker and with every loop I change the t_range and the dt. I observed that the memory is increase with every loop. I wanted to find out why this happened, so I broke the code down to the bare minimum. Therefore, the example above is only meant to show that the memory is increased with every iteration imitating the structure of the code without having a deeper meaning

[david-zwicker]

I see. I'm not surprised that memory consumption is increasing, since numba compiles quite a number of functions each time .solve is called. To avoid this, you could give the exact points where the tracker should hit to the interrupts argument, which would create a FixedInterrupts class. I don't know why you need to adjust the time steps, and suggest to use adaptive time stepping instead. If you need more control, you might need to write your own solver, using the ExplicitSolver class as an example.

[asmahr]

Thank you. But why does the memory not decrease again especially after using gc.collect()?

I will try the option to change the tracker. Currently, I am using the argument
tracker = [save_data.tracker(tracker_interval)] in .solvewhere, before the loop, I define the tracker_interval and the variable save_data = MemoryStorage()

[david-zwicker]

The point of caching is that the information is retained in memory. You might be able to convince numba to cache less, but this might be difficult to control in detail.

[asmahr]

I see

Since I am using already tracker = [save_data.tracker(tracker_interval)] and the memory increases, the only option is writing my own solver with backend = numpy to reduce the memory as you have proposed or is there anything else you suggest me to try?