Hey everyone,
Thank you for all of the kind words!
You can check out what my first draft of the end game (finite elements for unstructured grids to solve stochastic partial differential equations) looks like at this tutorial page:
https://juliadiffeq.github.io/DiffEqDocs.jl/latest/tutorials/femstochastic_example.html
It’s completely usable and you can already do some cool stuff, but let me say up front that this part is only a first draft and not well-optimized for finite element PDE solvers (I mentioned that in the blog post as well). In fact, I wrote this before most of the ODE solvers, and so it served as a jumping off point for how to “do this better for a 2.0”, namely:
- It needs to be changed to inplace updates instead of the vectorization I used.
- It needs to utilize ODE solvers. I actually wrote each method directly in here. The reason was because not only could no ODE solver do splitting, but no ODE solver could use the fact that the Laplacian is a linear operator, and thus when split we need to solve a linear problem. The naive splitting for Crank-Nicholson (Trapezoid rule) in ODEs actually requires a nonlinear solve.
- It needs to have a way for users to choose how to solve the nonlinear and linear problems. I had to implement a bunch of options to make this pretty extensive.
- Choosing algorithms by symbols results in dynamic dispatch. It’s not a big deal since most of the time is spent in an inner function and so the dispatch cost is fairly minimal, but it’s still there.
If you notice, all four of these problems were addressed by the design of the ODE/SDE/DAE solvers, and so when the setup with all of the special problem types and promotion is complete, the finite element solvers will switch over to just discretizing the problem, creating the matrices for the operators, and calling the ODE solvers (giving users the option for which algorithm). Since they will be able to recognize linearity, then for example Trapezoid() would directly give the Crank-Nicholson method with linear solves you know and love for the Heat Equation and all of that jazz. In fact, the event framework for ODEs/SDEs/DAEs is general enough that you can use it to implement adaptive meshing as well: use whatever criteria you want to trigger a meshing change, and the events can change the size of the ODE and replace the linear operator. This case will need a slightly different solution object to save a timeseries of meshes as well, but that’s not difficult.
But this is why the API for finite element methods is not quite up to speed with the other solvers: it’s old and there’s no reason to update it until this big change occurs.
All of that aside, a finite element problem discretizes into an ODE on matrices, where each column vector is the spatial discretization of a variable. So sol[i] is the matrix of values at the ith timestep, sol[i,j,:] are all of the variables at the jth location in space, and you have sol.mesh.elem[j] which gives [n1 n2 n3] as the nodes for the triangles, and sol.mesh.node[n] gives the (x,y) (and in the future z) coordinates for the node n.
So all of the information is there, but there’s probably a little bit of work to be done to make it a little bit more accessible. One thing which goes a long way here is plot recipes: the plot recipes take this and automatically make trisurf plots and animations. Note that there is a problem with Plots.jl that its trisurf doesn’t let you choose what node triples create triangles (like PyPlot does with its triangles argument), so the plots don’t look so great right now, but that issue was already noted a long time ago and fixes will come.
But together, this package shows a way to do finite element discretizations, and I’ll use it to provide some linear finite element discretizations. That should also serve as a good guideline for how other finite element discretizations could be done, which would have their own problem/solution types acting in a similar manner.
