Stream: helpdesk (published)

Topic: Advice on DAE-like heat transfer problem, DiffEq.jl

waterlubber (Sep 02 2026 at 01:36):

Summary: Should I model a DAE as a proper implicit DAE, use a "lagged" copy of the previous dX/dt term, or restructure my code to have much more state, but a simpler structure?

I am looking to simulate a "packed bed" heat exchanger, using OrdinaryDifferentialEquations.jl. I've used this a few times in the past, but have limited background in advanced math (just a 4-yr MechE degree), so

My model is currently a 1D ODE in time; a flow of gas passes through a packed bed, and I intend to simulate the thermal transient as it exchanges heat with the bed.

Currently, my state for the system contains the temperature of the gas and the bed packing. Calculating the change in gas & bed temperature over time is easy; however, this is dependent on the velocity through the bed, which varies according to density changes.

Unfortunately, this velocity term is dependent on the rate of change of density; and thus temperature; of the bed; my state vector dX/dt is a function of both X and Xdot. Thankfully, this effect is fairly small (the process is heavily advection dominated, so the effects of density change from temperature are likely on the order of ~5-10% of the total mass flow).

I am considering the following solutions:

1. Model the system as a proper implicit DAE by augmenting the state vector. From what I've read, DAE solvers are much slower, much stiffer, and generally look to be WAY over my head in terms of theory, so I'm hesitant to reach for this. However, if having a "small amount" of algebraicness (?) on what is otherwise a plain ODE is still performant, I am happy to consider this option. I do not see an easy way to state this problem in mass matrix form; all of these terms are very nonlinear and interact in tricky ways.
2. Simply iterate the solution in the solver until it converges.
3. Use a Jacobi-like approach; grab the previous state, pull its dX/dt, and use that to calculate the correction term. I am confident this will work, since I've solved very similar problems to this one with extremely simple first-order-Euler type Excel solvers. However, DiffEq.jl uses much more advanced solvers and I would need to be very careful about re-using the previous state. I think it's doable with a callback function, but I'd also need to force the solver to take small timesteps to keep this valid.
4. Re-write the code to hold mass as part of the state vector; instead of holding pressure constant and allowing mass flow to propagate instantly, treat the system as a series of connected nodes and allow the system to come to equilibrium in the way of any normal CFD solver. This seems like it would be slower, but it's conceptually simple and might avoid the severe penalties of DAE.

I am very new to advanced numerical methods here, and Julia as a whole. I've written most of my previous solvers in Excel with really basic Euler (or occasionally RK) integration.

Mostly looking for practical advice, or just additional context I might be missing. I'd love to learn more about the subject rather than just brute forcing a goofy solution.

(If this is posted in the wrong area, I sincerely apologize. More of an application question than an issue with the DiffEq library itself).

Daniel Wennberg (Sep 02 2026 at 12:31):

The structure of your equation is not entirely clear to me. It would be helpful if you could show a simplified example of an equation with the same structural issues.

Is your equation actually fully implicit, i.e., can $\dot{\boldsymbol{u}} = \boldsymbol{f}(\boldsymbol{u}, \dot{\boldsymbol{u}})$ only be solved iteratively? Even if it looks implicit when written like that, it may still be possible to solve explicitly for $\dot{\boldsymbol{u}}$ using forward substitution, i.e., there may be a "triangular" structure in $\boldsymbol{f}$'s dependence on $\dot{\boldsymbol{u}}$:

$$\begin{aligned} \dot{u}_1 &= f_1(\boldsymbol{u}) \, , \\ \dot{u}_2 &= f_2(\boldsymbol{u}, \dot{u}_1) \, , \\ \dot{u}_3 &= f_3(\boldsymbol{u}, \dot{u}_1, \dot{u}_2) \, , \\ \ldots \end{aligned}$$

such that you can eliminate $\dot{\boldsymbol{u}}$ from the right-hand side:

$$\begin{aligned} \dot{u}_1 &= f_1(\boldsymbol{u}) \, , \\ \dot{u}_2 &= f_2(\boldsymbol{u}, f_1(\boldsymbol{u})) \, , \\ \dot{u}_3 &= f_3(\boldsymbol{u}, f_1(\boldsymbol{u}), f_2(\boldsymbol{u}, f_1(\boldsymbol{u})) ) \, , \\ \ldots \end{aligned}$$

If this is the case, you can just treat it as any old ODE.

Apologies if you've already concluded this is not possible. It just seemed from your description that this might be the case, i.e., that you might be able to do forward substitution from $\dot{T}$ (temperature) to $\dot{\rho}$ (density) to $\dot{\boldsymbol{v}}$ (velocity) or something similar.

waterlubber (Sep 03 2026 at 01:04):

Thank you Daniel.

I believe the equation is fully implicit.

To better summarize my system:

The system is a series of cells; one set represents fixed packed bed material and another set the gas.
Each timestep, heat is exchanged with the packed bed material (classic ODE style), and advected along the cells.
Mass flow out of the rightmost cell is held constant (set by other simulation parameters, actually); this is then propagated left using the density change of each cell to compute the next mass flow.

However, the density change in each cell is dependent on the previous state derivative dT/dt; this advection term is pretty significant at the boundary and thus the entire simulation is dependent on the previous derivative. Furthermore, the heat transfer term _between_ the bed and gas is dependent on mass flow, although this dependence is much less than the effects of advection.

What you are describing with the triangular structure looks a lot like what I called "lagging," pulling the previous state's dT/dt value and using it to compute the mass flow values for the current step. On a whim, I set up the simulation with a first-order Euler approach and very small timestep, exactly as you described. However, it still exhibited instability; when hot gas flows into the first cell, the advection term (warming effect of hot gas) produces a very large, but real, dT/dt value. The corresponding mass flow from this rapid change in density is negative and causes the solver to oscillate. I don't think there's a way to solve this without some sort of implicit solver, where the input mass flow was reduced until the system converged.

That leaves me with a few remaining approaches:

1. Using a proper implicit DAE solver
2. Adding a mass state to the gas cells, and allowing pressure to float. Mass flow would then be computed from the difference in pressure between cells. This difference would be very small (< 1% of the cell magnitude) and I suspect this approach would be fairly inefficient.
3. Implicitly solving for mass flow within the derivative solver; i.e iterating the mass flow into the first cell until the process converged. I don't know if this would be more efficient than a plain DAE solver.

I plan to add a lot more scope to this simulation; I think this is the only part of the entire program that would require a DAE.

Appreciate any advice you have to offer, and thanks for responding in the first place.

Daniel Wennberg (Sep 03 2026 at 15:42):

What exactly do you mean by "previous state derivative dT/dt"? Is there a finite time delay here? (In which case you have a delay differential equation, DDE.) Or is this just an artifact of intuitively building up the model in a discretized fashion, and once you take the continuous-time limit you obtain the instantaneous dT/dt?

waterlubber said:

What you are describing with the triangular structure looks a lot like what I called "lagging,"

There's no lag/dependence on previous derivatives in what I wrote. The indices on $\dot{u}_i$ denote different components of the state vector $\boldsymbol{u}$, not time steps or anything like that. I was just describing a case where it's straightforward to take an ODE that looks like it's implicitly defined ($\dot{\boldsymbol{u}} = f(\boldsymbol{u}, \dot{\boldsymbol{u}})$) and solve for $\dot{\boldsymbol{u}}$ to make it explicit ($\dot{\boldsymbol{u}} = g(\boldsymbol{u})$).

I do feel like something like this should be possible for your system, i.e., solving explicitly for all the derivatives from right to left, if you include enough state in your state vector. Something like the mass or the density in each cell. Perhaps you'd be left with one algebraic constraint for conservation of mass, but if I understand correctly, that's an index 1 DAE on mass matrix form, which is much less intimidating than a fully implicit ODE. (Caveat: I'm not at all an expert on DAEs and such, I'm just a fairly experienced user of DifferentialEquations.jl and peruser of its documentation.) You say you suspect such an approach would be fairly inefficient, but it may be worth at least giving it a try. An explicit equation in a larger dimension may be preferable to an implicit equation in a smaller dimension, both in terms of raw efficiency and because it gives you access to a much larger array of solvers with more features.

If you end up sticking to a fully implicit formulation, I would definitely go for a fully implicit DAE solver (Sundials.IDA seems to be the go-to) rather than setting up an iterative solve within your ODEFunction. The latter combined with adaptive timestepping sounds like a terrible idea for efficiency, as you'd be solving $F(\dot{u}, u) = 0$ from scratch over and over and over. A DAE solver combines timestepping and the $F(\dot{u}, u) = 0$ solve in a single system of equations and can reuse Jacobians across timesteps,.

waterlubber (Sep 03 2026 at 18:06):

It is an artifact of me building the model in a discrete fashion; I have very little formal education in the way of differential equations and most of my exposure through the field was through a course on numerical methods and CFD. Apologies for this; I mostly end up thinking of problems like this with real-world values and discretizations of everything.

I can actually rewrite my system in such a form to permit this; this is by including a density term for each cell. This lets me remove the direct algebraic constraint for constant pressure by allowing the pressure & mass in each cell to "float"; essentially it's spreading out the iteration required to solve over multiple steps.

I have set up exactly this right now and am experiencing the expected stiffness/stability issues (namely, the mass flow between cells is based on the pressure drop between cells; this is very small compared to the magnitude of the pressure and there is catastrophic cancellation).

I might try a fully implicit method on the weekend to see how it goes. In the meantime, I might just hack this part of the system out and replace it with an unphysical model that replicates the qualitative performance.

On the topic of DAE solvers: the vast majority of this system is a plain ODE, with a very small algebraic component. I do see that the DAE solvers require an array that indicates differential components for initialization.

Will the solver be able to take advantage of this structure? (i.e, small algebraic component largely uncoupled from the rest of the system?) Is there any hinting I could do that might make the process more efficient?

Update: the explicit, pressure based method with extra state oscillates like crazy. Probably would need to use an implicit solver for this.

waterlubber (Sep 03 2026 at 19:02):

As an addendum, a half-remembered technique from my CFD class for similar unstable advection dominated problems was leapfrog integration. I'm not sure if it's directly applicable to this case but might provide a useful starting point for others with similar problems.

Daniel Wennberg (Sep 03 2026 at 19:13):

I would assume that DAE solvers are designed to do the best they can with the structure you give them, though I don't know the details of how they work internally.

Not sure I have any more insight to contribute, but I hope you can figure out a solution. Let us know what you learn!

waterlubber (Sep 04 2026 at 18:11):

I ran with the Sundials IDA solver and after fixing all my variable names / refactoring overhead, simulation converged on a realistic looking result on the first try. Performance isn't great, but it's still acceptable - takes about 10 minutes to run with 64 cells for a minute of real runtime. I imagine most of the overhead is in the super slow PropsSI() from CoolProp, and I might be able to improve it by swapping to Clayperon.

Thanks for your help! Leaving this up for any future people in a similar boat. The Sundials IDA solver is Just Really Good and should be your first try if you have freaky implicit equations.

Daniel Wennberg (Sep 04 2026 at 18:36):

One thing I wanted to mention, have you tried using ModelingToolkit.jl to build your model declaratively and let it handle all the lower-level stuff? I haven't used it much myself, but it's supposed to have a lot of sophisticated functionality for things like structurally simplifying DAEs. Maybe not so well suited for a large, finely discretized model (unless you just give it the underlying PDEs and let it work from there), but would be interesting to try it on a much coarser discretization just to see how it handles the structure.

waterlubber (Sep 04 2026 at 18:48):

I took a cursory look at ModelingToolkit, but I suspect PropsSI would break it entirely. CoolProp calls out to a bunch of non-Julia code and all the derivatives are computed numerically (it barely has Unitful support, let alone autodiff) so I don't think I'll get a lot of symbolic simplification.

Christopher Rackauckas (Sep 12 2026 at 12:26):

For such a thing, DNordseickBDF should be faster these days

Christopher Rackauckas (Sep 12 2026 at 12:26):

yeah the coolprop part can usually be slow

Christopher Rackauckas (Sep 12 2026 at 12:26):

did you profile to see where your bottleneck is?

Last updated: Sep 19 2026 at 08:53 UTC