Time-dependent problems with Method of Lines using JuliaDiffEq

I’m writing a code to integrate a time-dependent PDE system that I convert to a large ODE system using Method of Lines. These are very long time integrations, for now in 2+1 dimensions, which I’ll eventually generalize to 3+1 problems.

I have used JuliaDiffEq for simpler problems (in 1+1 dimensions) and it works great. However, for my present task in which I’m dealing with very long simulations that can run for days or weeks, it seems that JuliaDiffEq is not very suitable. Is this the case? If JuliaDiffEq is actually suitable for this task, is there any example that I could have a look at?

You optimized your code and integrator choice? I would run the profiler on your derivative function and make sure it’s all good, and test out a few integrator methods. If you need a stiff integrator on a large system, make sure you aren’t using a dense Jacobian unless it’s actually dense.

I have a simple RK4 integrator that I was planning to use for this, once I have everything in place. I’d be happy to use JuliaDiffEq if I can, though. The thing is that puts me off is that it seems that JuliaDiffEq is designed in a way where one sets up the problem which is then solved in one go through sol = solve(prob).

This won’t work for me since there is no way my whole integration will fit in memory, so I’ll need to write the output as the integration moves along and forget about the previous timesteps. Is there a way of using JuliaDiffEq in this fashion?

See the integrator interface. I assume what you’re looking to do is query to step as you please and then serialize outputs at specific times? That will do it. You could also use the FunctionCallingCallback

Ah, this seems nice, thanks. Yes, I want to step as I please and output at specific times but also, crucially for memory constraints, as the integration marches along I can only keep data from the current and previous time step. Otherwise I’d run out of memory pretty fast… Can this be done as well?

Yes, save_everystep=false. You may want to go through the common interface page which explains the saving arguments. It saves as much or as little as you want.

I’m not on a plane anymore so it’s easier to link. This part of the common solver options page describes the output controls:


Here’s some quick examples:


You can also completely take over saving by turning it all off and using the integrator interface


or by turning saving off and using one of the callbacks:


Continuous dense output is the default because when people solve most simple ODEs that’s the nicest output, and that tends to be the majority of users first learning DifferentialEquations.jl. However, as you get to harder problems, you usually want to make a conscious choice of how you’re saving. That’s impossible to automate so we just make sure all of these options are there.

perfect, thanks!