JacquesCarette / Drasil

Generate all the things (focusing on research software)
https://jacquescarette.github.io/Drasil
BSD 2-Clause "Simplified" License
142 stars 26 forks source link

Infinite recursion possible in equations exposed by `DifferentialModel`s formed for type-checking #3322

Open balacij opened 1 year ago

balacij commented 1 year ago

image

If we have a 1x1 matrix of coefficients and only 1 unknown, then the formed equations for type-checking can be an infinite list, resulting in an Out-Of-Memory (OOM) error.

A highlighting code prototype is public. Run make bmbnd_diff to reproduce from that tree.

gotop:

image

smiths commented 1 year ago

@balacij I can't easily see what ODE you are trying to create. Is it possible that the ODE is nonsensical? (This should still be detected of course) I don't really know what you are trying to do in Drasil with this example.

balacij commented 1 year ago

It may very well have been nonsensical. I was trying to encode (E_B * I_B) * d^4 y_B/dx^4 = w_B(x). On second look, w_B(x) is definitely not a 'constant'. Either way, as you mention, the error should be detected if it is.

https://github.com/JacquesCarette/Drasil/blob/b73f49ff4a144322d1b1840d43a45b4ad64a074e/code/drasil-lang/lib/Language/Drasil/Chunk/DifferentialModel.hs#L160-L177

smiths commented 1 year ago

The ODE makes sense as written. Would it help if we manually converted the fourth order ODE into a system of 4 first order ODEs?

balacij commented 1 year ago

That is an interesting idea. After a bit of exploration work with the existing used ODE solvers, I think it would only work if we can also "erase" the boundary conditions. The existing ODE solvers are all focused on providing initial values for the dependent variable at an initial value for the independent variable. Here, since we're looking to encode something with more than boundary values as well, we would not have that included/considered in the solver. Do you know if we can somehow include the boundary conditions in the system of equations to circumvent this restriction?

smiths commented 1 year ago

No, I don't think there is a way to include the boundary conditions in the system of equations. The boundary conditions are needed to find a unique solution, but they are handled differently than the system of ODEs.

We should rely on a library to do the hard work for us. We just need to get the information to the library. The BVP solver you found from scipy should be fine. We need to tell the BVP solver the ODE and the boundary conditions.

I don't think you need to worry about a system of ODEs. My understanding from talking to you is that the scipy solver takes the fourth order ODE directly. If we do need to turn it into a system of ODEs, we would have something like:

EI dc/dx = w(x) a = dy/dx b = da/dx c = db/dx y(0) = 0 y(L) = 0 b(0) = 0 b(L) = 0

(The last two boundary conditions are from d2y/dx2(0) = d2y/dx2(L) = 0.)

JacquesCarette commented 1 year ago

This will require extending ModelKind for the ODE case to signal that it's a BVP. These are different enough that they deserve their own case. That will then help drive the back-end's choice of a BVP solver.