DifferentialEquations.jl FEM example

General Usage
1

I want to solve a time-dependent PDE on a somewhat nontrivial geometry. As a step in this direction I am looking for an example how to solve a PDE using DifferentialEquations.jl + some FEM package. For instance, the heat equation or a simple advection would be interesting.

2 Likes
2

We used to have examples like this, but they were never updated. Basically, just use FEM to get an ODE and then solve the ODE. Most packages like FEniCS (FEniCS.jl) will kick out the matrices as sparse matrices to build said ODE.

2 Likes
3

Thanks! Can you maybe point me to such an example, even if it won’t run anymore? Here is what confuses me:
For DifferentialEquations.jl I need to encode the problem as problem(du, u, p, t). For his to make sense du and u must live in the same vector space. Now it seems to me that FEM spaces are not stable under e.g. taking derivatives. (Variational formulation is used to circumvent that problem). So it is not clear to me how to use FEM to bring the problem into a form consumable by DifferentialEquations.jl.

4

FEM turns problems into ODEs for the coefficients (at least for time dependent discretizations), possibly with a mass matrix. So write it in Mu’ = f form. du and u technically never live in the same space: they have different units in all physical problems, but that’s not an issue and is just a part of the definition of f.

1 Like
5

Thanks. That ∂ₜu and u have different units is a good point. However I still don’t see how to get the PDE in problem(du, u, p, t) form. Can you maybe explain what the Mu' = f form is in case of the heat equation ∂ₜu = Δu?

6

For the heat equation, you generate a stiffness matrix A and a mass matrix M when you do FEM, and then you just do Mu' = Au.

7

FinEtoolsHeatDiff.jl implements a simple heat diffusion transient problem in the file hill_decay_diffeq.jl. The code runs, but the results are weird: there are some initial oscillations, which shouldn’t be there, but the temperature results towards the end of the time span are roughly correct (@ChrisRackauckas is hopefully looking into it). The file hill_decay_gentrap.jl implements the hand-coded implicit Euler solution of the same problem.

To run, activate and instantiate the package and include the files.

EDIT: Found an error in my code. (DifferentialEquations.jl was fine, all my fault.) The result is now essentially the same for the hand-coded version and the DifferentialEquations.jl form. One more optimization is possible: using the mass-matrix form of the solver. That doesn’t work at the moment, and we (Chris and I) don’t know why.

EDIT 2: The mass-matrix form now works as well (https://github.com/PetrKryslUCSD/FinEtoolsHeatDiff.jl/blob/master/examples/transient/2-d/hill_decay_diffeq_mass.jl).

1 Like
8

your problem(du,u,p,t) has to assign the spacial derivatives to du:

let’s assume for simplicity a equidistant 1D-mesh (the mesh is stored in the parameter p. Your function has to return one equation per stencil ( telling the solver the ∂ₜu at each grid point)

function problem(du,u,p,t)
 du[1] = 0   #Bondery condition1
 for i in p.stancils[2:end-1]
   du[i] =  ( u[i+1]-2*u[i]+ u[i-1] ) / Δ^2 #second derivative of u
 end
 du [end] = 0 #Bondery condition2
end

This function then is your f in Mu' = f , and the du is the Mu'. So whats so special about the mass matrix M? In case of fixed Dirichlet boundary conditions it is possible to define the problem solely by the indication of the temporal derivatives. In this case M = I (which is the default value when no mass matrix is specified). But if you look at the function above it is possible to multiply the whole PDE with Δ^2 to get rid of the division, resulting in a mass matrix of M= I * Δ^2 . Finally you can express the function problem(du,u,p,t) as a matrix matrix multiplication of

 A = 
[0  0  . . . 0
 1 -2  1 . . . 0
 0  1 -2 1 . . . 0
.  .  .
0  . . . 1 -2  1
0  . . . 0  0  0]

function problem(du,u,p,t) 
 du = A * u

and write Mu' = Au

How ever this is just an example which is far from being performant.

9

I’ll look tomorrow.

10

If it is just heat diffusion you could try mass lumping to get rid of the mass matrix.

11

I’m possibly planning to try and implement an existing FEM model that I have using JuAFEM into DiffEq, if I finish it I’ll post the result here as a possible example.

12

True. The price to pay is the loss of symmetry of the “Jacobian” (df/du).

13

It’s not difficult though:

ff = ODEFunction(f,mass_matrix=M)
ODEProblem(ff,u0,tspan)

etc. Just a bad day to respond haha.

14

Fact is solving time-dep PDEs is usually best done with taylored hand-coded algorithms. General-purpose packages are of course fine for prototyping. Performance is usually a different matter.

15

Completely false. In fact, most literature shows exactly the opposite of this. It’s very small systems, like small SDEs, where you can eek out advantages, maybe a few linear cases (though most differential equation solvers have options to tell it that it’s the linear case). For PDEs, the op overhead is so high that pretty much none of the small details matter (one matmul costs more than all of the scalar exponentials of an adaptive solve combined!), and what really matters is using a strategy that minimizes the number of right hand side evaluations. Try to “optimize” some second order time stepping routine is pretty much the biggest waste of time, since OrdinaryDiffEq, Sundials, Hairer’s routines, etc. will already be non-allocating during the solve, and you can even see that at the size of PDEs allocations basically don’t matter (again, O(n^3) matmuls are the same reason why ML just uses functional programming styles).

What makes PDE solving fast is using things like Numerov schemes, Rosenbrock-W methods, adaptive time adaptive order BDF, IMEX RK, ExpRK, etc. with additions like PI-adaptive time stepping and tuned Jacobian reuse strategies. These methods take less time steps, less Jacobian evaluations, and less matrix factorizations than any of the simple Crank-Nicholson or whatever implementations out there. This is why all of the main high performance PDE solving libraries like PETSc have spent so much time tuning the stepping strategies, and it’s why optimized HPC PDE solvers are not using handwritten time stepping codes (with the most common calls being to Sundials and PETSc of course).

(And probably the weirdest thing is when people point to “PDE specific” methods. I think the biggest shame in the pedagogy of PDEs is actually how little it is taught that the “PDE schemes” are actually methods on specialized ODE formulations. Exponential integrators are just methods on semilinear ODEs. IMEX RK methods are just on f = f1 + f2. Even ADI is just an approximate matrix factorization on Trapezoid. This insight makes it easy to develop a method which dominates over “standard” hand written methods, like do approximate matrix factorization (“ADI”) but in the context of the 4th order IMEX KenCarp4 with PI-adaptivity and with a good reuse heuristic. See how ADI compares to that…)

4 Likes
16

Obviously our experience differs. But Chris is welcome to put some numbers behind this theorizing.
https://github.com/PetrKryslUCSD/FinEtoolsHeatDiff.jl/tree/master/examples/transient/2-d
If you can, run the sim faster than the hand-written version.

17

Where is the work-precision diagram? You’re solving with OrdinaryDiffEq’s methods at 1e-4 tolerance whereas your implicit Euler is solving with dt=2.0, which we know from the previous experiments is at a lot higher error than 1e-4. If you actually compare method timings to get to the same error, not a simple method with high dt vs a higher order method getting a much lower error, and do the obvious tweaks (you’re using Cholesky, why not setup the linsolve to Cholesky instead of lu?) you see a major difference. There is no theorizing here, this has been shown 100 times already…

(What’s going on here is you specialized on the linearity. Make it nonlinear, or set the same specialization in DiffEq)

18

It is easy to compare by showing absolute error vs. solution time. Again, why not show actual numbers?
I haven’t done any optimizing, but the hand-coded version runs in around 0.01 second and the mass-matrix form from DifferentialEquations.jl runs in around 3.3 seconds. If we go through a proper assessment of accuracy versus solution time, then we will know more. But these are informative numbers.

EDIT: As Chris pointed out, these numbers should not be taken too seriously as the DifferentialEquations.jl solver runs with a dense Jacobian. It also apparently does not recognize linearity of the differential equation. I am currently trying to understand how to incorporate these features so that the comparison is fair. The work-precision diagram is one way. But the cost needs to be calculated as they say comparing apples to apples.

19

Yes, so make a WP plot of the time series against a reference solution computed at a very low tolerance.

I told you how to do it… I’ll post it when I find the time to do it here. I don’t know why this is so hard to understand, but since this has been done so many times it’s not really a priority. I may find time over the weekend.

20

The issue in your example is that your hardcoded implicit Euler is using a sparse Jacobian while you let DiffEq default to a dense Jacobian. I don’t see why you would think that a dense 2601x2601 would be as fast as a sparse one? Why not… just set the options?