New user here, quite eager to get involved with the community, so my apologies for the long post!
I work in fusion research, and am interested in testing modelling approaches for systems of coupled stiff integro-differential equations (with stiff ODEs and PDEs in the mix). The equations in question are kinetic and fluid equations coupled with Maxwell/Poisson equations and collisional-radiative models in general. I’m new to Julia (coming from Modern Fortran), and have been looking into existing packages (DifferentialEquations, ModellingToolkit, etc.) to figure out what sort of methods are available. Coming from a low-level language, I’m quite happy to see a lot of easy to use packages. I would ideally like to have a framework where I could have control over both the equations I’m solving, as well as the discretization details.
My (infinitely ambitious) wishlist for this sort of framework would be something like:
Be able to select different discretization options for different spaces in the domain (e.g. discretize the spatial coordinate using one method and the velocity coordinate using another)
Have control over the time-integration so that I could have different variables evolved using different methods (implicit vs explicit) as well as have Strang-like splitting of terms in some of the equations (this is important for the stability of some kinetic equation discretizations).
Be able to work with highly non-linear stiff systems
Have a choice of various discretization options (finite difference, volume, element, as well as particle methods)
Not shooting myself in the foot in terms of performance
The naive Modern Fortran approach would be something like:
Write complicated variable objects to handle the requirement that different variables might be defined on different grids/spaces
Suffer through MPI to make sure that halo-exchange/communication works at an early stage
Construct abstract operator objects that can transform my variables and combine them into equations
Build a time integration system that can handle the above requirements (likely using PETSc or something similar)
Pray to whatever deity I can think of that I don’t destroy performance through abstraction
While I’m sure I could go through the above naive approach using Julia, I also feel that it would be unnecessarily painful and unproductive (though likely not as painful as Fortran). Looking at the current state of ModellingToolkit, a lot of these features seem to be on people’s minds, though I’m unsure of the exact state (I have the JuliaConn ModellingToolkit Workshop on may watch list).
I would appreciate any feedback from the community, especially in terms of existing tools. I also realize that my ideal framework is likely impossible given its many requirements, but I wanted put it out there nonetheless.
For ordinary integro-differential equations of specific forms, we may want to develop specific problem types (IDEProblem etc.) so that we can write specific solvers to the specific forms. We just haven’t gotten there yet, but it’s part of the DifferentialEquations.jl general plan in some sense.
Thanks for the reply! Great to know that this is on your radar.
My main interest is in exploring numerical methods (and combinations of numerical methods) for solving kinetic equations, so I guess that I would have to build my own Julia type for the variables and the specific sparse matrices for different discretizations of the collision operator (for things like the Boltzmann equation). Given the specific goal of trying out new numerical methods I don’t think I can start at a level higher that this, correct?
I’m looking forward to exploring DifferentialEquations.jl for solving the resulting ODE system. I did see there is basic linear-nonlinear splitting for IMEX methods, but I couldn’t find anything more general that would let me solve different equations in the system using different methods while keeping them coupled. I might not have looked hard enough, though.
I think “type per PDE” will lead to non-composability, hence the MTK push to having a single PDESystem with symbolic representations. However, in the IDE space I think there is a general set over which methods can be written, so we just need to formalize that set itself and look at how other methods are specializations of one general form.
Ah, I think I didn’t express myself properly. What I meant was not to build a different type for each problem (PDE), but to build a struct that I could use as an initial condition for ODEProblem, and then write my own update functions that would use operators based on some sort of sparse matrix representation. I think I saw another post somewhere around here about using user-defined types as initial conditions for ODEs, so I wanted to try that, but with me doing the matrix building. Is this still an ok approach or is there something in MTK that makes it obsolete?
I probably should have mentioned that the problems I’m interested in have non-trivial BCs and likely complex domains/geometry, so I’m not sure whether I can avoid having to manually build at least some of the sparse matrices, given that I’d like to have access to different discretization approaches.
I’ve just seen that you can compose symbolic equations in MTK. Is there any plan to do something like this for solvers, i.e. choosing a different ODE solver for each of the equation subsets?
Thanks for your patience, I’m still wrapping my head around the Julian way of thinking!
Wrt choosing different ode solvers per equation group, I think that’s basically the goal that IMEX solvers are working toward. Having more flexibility here would be great though, and if you want to help with that, one simple thing that would be good is getting sample problems that show a benefit from this approach.
The main class of problems I was looking at as candidates for the solver splitting are coupled advective-diffusive-reactive equations where only a subset is stiff. While I haven’t tested this, I suspect that a partially ionized plasma would be such a case. Another plasma example is solving the Vlasov-Maxwell system with the fields evolved implicitly and the plasma explicitly - allowing you to skip resolving the shortest timescales. I’m not sure either of these problems is simple, however…
Another solver feature that I’m struggling to find in the DifferentialEquations documentation is support for Strang splitting. Is this (or any other time-splitting method) implemented anywhere?