Julia's FEM Landscape: where should DifferentialEquations go?

Specific Domains Numerics
1

Hey, I am taking a good hard look at the Julia FEM landscape and would like some input. Let me first explain what I see the domain and scope of DifferentialEquations.jl. My work and research is centered around solving differential equations, but mostly related to timestepping routines (for both continuous and discrete stochastic equations). I plan to stick to timestepping routines because I am making a lot of progress and nobody else seems to be doing it (and it’s more interesting of course :slight_smile: ).

So while this involves PDEs (/RDMEs), there is a lot of work that has to be done in meshing, assembling, etc. that I would like to make use of instead of re-implementing. But I see DifferentialEquations.jl as continuing to offer:

  • Libraries of timestepping routines for ODEs/SDEs/DDEs/DAEs
  • Well-optimized end-to-end solvers for common PDEs.

While before I internally created a small FEM toolbox to offer part 2, continuing in that direction seems impractical. Instead, I am looking around to see what’s the right FEM toolbox, and working with those who are making this offering. That said, I don’t know what’s the right choice here.

On one end, I think that plugging into some established FEM suites like deal.ii and FENICS might be a good idea, but there seems to be no work done on developing interfaces to these in Julia? I would be weary of such a solution because I wouldn’t want to lose the type-genericity though: since the most optimized timestepping routines are the native Julia methods, the full set of Julia-defined numeric types are supported, and so if the FEM toolbox doesn’t support the full type system (and is restricted to Float64) it would be a big loss.

But I am not sure about the native Julia offerings. @kristoffer.carlsson 's JuAFEM.jl looks to be an interesting choice (which offers assembly), but I am not sure what the long term plans are there for integrations with meshing tools the different types of elements you plan to offer (mesh adaptivity as well?). On the other hand, I am not sure what JuliaFEM offers. It seems to be more along the lines of FEA tools which are built for physical applications rather than a toolbox of primitive functions for a standard FEM workflow? I could be mistaken.

If anyone knows this area and Julia’s developments in this area well, I could use the guidance. Thanks in advance.

11 Likes
2

This is pretty good statement of the current status. However, we planned to split the package to smaller pieces for v0.3.0. At least all io related stuff will go out. Also we planned to promote the JuliaFEM organisation as a home for all Julia related FEM development. We are also open of the idea of renaming JuliaFEM.jl package, especially if it will be splitted (maybe meta for the JuliaFEM organisation).

Anyway our original scope was reliable modular tool for real industrial problems is still valid for us. We haven’t reach the version v1.0 yet but still we have something working which have created value.

Can you define what you mean by this? Any example even from another open source FEM is enough. In my experience this means different things for different education background.

Few academic examples:

3

I’m thinking about something like FENICS or deal.ii: they provide tools that bring you through meshing and assembly and leave you with this discretized version of the differential operators.

This is a good example from FENICS:

It has all of the tools to define the equation in terms of discretized differential operators. I want to have the backend do just this: use an FEM toolbox to build up some standard discretizations (these steps are shown as lines 1-46), and then solve these problems with various timestepping routines (for time-dependent problems, pluging in some linear/nonlinear solvers as needed).

4

I’m not sure what your plan is. You want to provide solvers for “common” PDEs but isn’t the simplest to just refer people to the FEM packages that already exists then?

Just some comments about JuliaFEM vs JuAFEM.

JuliaFEM and JuAFEM have different purposes. JuliaFEM seems to be more similar to a commercial full fledged FEM solver where most things are automated for you. JuAFEM is more similar to deal.ii where you implement your own weak form and do the assembly and solving yourself. JuAFEM just tries to simplify things like evaluating shape functions and nodal functions for different interpolation spaces for different quadratures etc.

As a fun side not I am currently taking a course from the original creator of FeniCS (Anders Logg). It is quite powerful when you have a standard weak form that you want to solve on the whole domain but if you get more complicated stuff like plasticity it seems that most of the nice abstractions break down.

1 Like
5

Yes, people always can just use FEM + ODE/DAE/DDE/SDE packages to solve whatever PDE they have, and I would never try to offer “every PDE” or something silly like that. But some disciplines spend a lot of time solving a very specific equation. In systems biology, various forms of semi-linear Heat equations (Reaction-Diffusion equations) are very common, with a lot them being stochastic/stiff/etc., and so there’s an extra step (an entire field of research) beyond “here’s a spatial discretization” which I develop methods for and might as well offer with an easy to use interface (insert reaction functions here, pick a mesh+BCs, and solve using timestepping algorithm(s) X).

I mean, I totally agree that this last step isn’t necessarily difficult to do implementation-wise, but the reason why I think it’s simple is because I’ve learned it and it’s what I do. A lot of people really do use things like the MATLAB PDE Toolboxes because it does lower the mathematical barrier to entry. And instead of just stopping 99% of the way there, I might as well offer an “insert function here” type of PDE solver, at least for some common problems which non-mathematicians and non-physicists are interested in (especially when I am already writing such software for my own research). I think this is adequate motivation to maintain such solvers (though they should require little to no maintenance, and would serve as a good set of integration tests anyways).

Okay, yeah, this is what I was starting to understand. JuAFEM is more like what I am looking for then. Do you plan on going all out on it? Like “feature-complete” with deal.ii/FENICS with adaptive meshing tools and all of that jazz, or do you plan to keep it simple?

1 Like
6

I would think wrapping FENICS would be a good first step. It has a large user base, who would be more likely to try Julia if they knew the FEM framework that they are already familiar with is supported. Support for types beyond Float64 is the kind of thing that you don’t realize you need until after you have it. (Even then, linear algebra beyond BlasFloat is still very rough in Julia.)

7

That’s why I’m wondering where JuAFEM/JuliaFEM is headed though. From what I’ve learned I think JuAFEM could be growing towards being a Julian drop in replacement/upgrade for FENICS. If that’s the case, I’d rather give it a try and help it along the way rather than building a wrapper that I’d be wanting to drop support for as soon as a Julian solution comes along. But if no one is working on such a package, then work on a FENICS wrapper makes sense.

8

A FEM package which would allow FeniCS-like high-level definitions but then also allow more low-level interaction (deal.ii-like) would be well cool. Julia would be very/uniquely suited for such a project. However, I suspect to get anywhere near feature parity with FeniCS/deal.ii will take a long time. So a FeniCS wrapper maybe a good intermediate step. (Just bike-shedding, I will not have time/expertise to contribute much to an effort like this.)

(Side note: a mesh library as proposed by Anders Logg (2012) and as is used in FeniCS was my first Julia project: https://bitbucket.org/maurow/lmesh.jl.)

4 Likes
9

This is a point that a lot of people have complained about. Having both approaches (Deal.ii, Fenics) within the same package seems very possible in Julia.

10

I think having some kind of domain-specific language in Julia for FEM assembly would be very nice. Our CFD code has such a system in C++ by using the Boost Proto library, but it is very difficult to maintain and to support students who have to develop new models. An example for the Poisson equation is here:
https://github.com/coolfluid/coolfluid3/blob/master/plugins/UFEM/src/UFEM/demo/PoissonProto.cpp#L56-L70

Or the stabilized incompressible Navier-Stokes equations:
https://github.com/coolfluid/coolfluid3/blob/master/plugins/UFEM/src/UFEM/NavierStokesAssembly.hpp#L46-L74

The idea is to assemble the linear system block by block, using epressions similar to the weak form of the problem. The advantage of this system is that the code becomes completely independent of the dimensionality and element type, and that the dimensions of the different element matrices are known at compile time, allowing in this case the Eigen matrix library to operate at full speed on these element matrices.

The main drawback of doing this in C++ template metaprogramming is that the code is hard to maintain. For example, I want to write an expression optimizer, but put it off so far because the task is too daunting. Also, supporting students who have to implement a model is not easy. To give an example: a compile error can result in a text file of several megabytes. I also measured the demangled string length of the longest type in the code, it was over 60000 characters long, so it’s no surprise that compilation can take over 4 GB RAM for certain files.

All of this is the main reason I’m moving to Julia now. I believe its powerful metaprogramming features make it the perfect to develop this kind of domain specific language, complete with sane error reporting and expression optimization and without sacrificing the C++ speed.

It’s a bit too much work to convert everything at once, so I intend to first wrap our current C++ code and then replace it piece by piece with Julia equivalents.

4 Likes
11

Yeah, this is why I think it has to be native Julia, because every time I mention the generic programming that I want, people mention that you can do it in C++ via templates but don’t do it because then your code is mangled and unmaintainable… so native Julia it is.

Do you plan on writing the Julia equivalent? It seems there’s more than a few people waiting on the sidelines for one to come along, but if no one is really stepping up to write it then a quick FENICS though PyCall wrapper sounds like an intermediate solution which we should consider.

12

Yes, but no promises regarding when :slight_smile:

From what I understand, FENiCS allows writing assembly expressions and then compiles these using a custom pipeline. I’m not sure if that makes it so easy to wrap, but if someone wants to try and it goes smoothly, it might be a quick way to get a mature package.

13

Using FENiCS with c++ you write the weak form in a DSL (basically python) in a separate file ke which you can compile to a header file ans then include in your project. If you use Python the weak form is compiled automatically when needed.

Wrapping the Python wrappers of FENiCS definitely seems like the easiest way forward.

14

I am working on a BoundaryElements.jl package. This package is centred around a limited number of fairly well defined concepts such as Space, ReferenceSpace, Operator, DiscreteOperator, … It offers functions for assembly of both integral and differential operators (so FEM assembly is a strict subset of the capability). For now the concepts have been implemented for the Laplace equation, Helmholtz equation, and Maxwell equations. General support also is included for transposing operators, systems of equations, the use of direct products of finite element spaces.

It is easy to extend the framework to include new operators and spaces (I have gone through this exercise myself whilst working on e.g. the 3D Helmholtz support).

For now, using Julia’s dispatch and macro capabilities I have designed a very rudimentary form compiler. For now it just provides syntactic sugar but the goal would be to build bespoke methods for assembly etc. JIT.

The goal and scope of this package are best demonstrated by a short example on how it can be used.

# Define the kernels
κ = 1.0
S = SingleLayerTrace(κ)
T = MWSingleLayer3D(κ)
I = Identity()

# Build the geometry
h1, h2 = 1/12, 1/15
h = max(h1,h2)
W, H = 1.0, 0.5
Γ1 = meshrectangle(W,H,h1)
Γ2 = translate(meshrectangle(W,H,h2), point(0,-H,0))
γ = meshsegment(W,W,3)

X1 = raviartthomas(Γ1, γ)
X2 = raviartthomas(Γ2, γ)
X = X1 × X2

# Define the incident field
E = PlaneWaveMW(z, y, κ)
e = (n×E)×n

j, = hilbert_space(:j)
k, = hilbert_space(:k)

trc = X->ntrace(X,γ)
dvg = divergence

α = 1/(im*κ)
β = α*log(abs(h))

Eq = @varform begin
    T[k,j] +
    α*S[trc(k), dvg(j)] +
    α*S.'[dvg(k), trc(j)] +
    β*I[trc(k), trc(j)] == e[k]
end

eq = @discretise Eq j∈X k∈X

b = rhs(eq)
Z = sysmatrix(eq)

u = Z \ b

fcr, geo = facecurrents(u, X)

All of it is written in Julia (I don’t rely on quadrature kernels written in c++, or existing functionality offered by e.g. FeniCS or BEM++).

As you can tell from the example, the package is light on the meshing and algebraic solving sides. Of course I will work towards (preconditioned) iterative methods (strengthened by matrix-vector accelerators such as the Adaptive Cross Approximation or Fast Multipole Method. There is descent support for time domain boundary integral equations. Actually Julia package WiltonInts84 contains tested and optimised routines for the fragile quadrature rules that form the core of these time domain methods.

I am planning to publish this and supporting packages in the following months if I can get them in a state that is acceptable. Very happy to talk/collaborate with interested people on this.

4 Likes
15

Sounds interesting. I don’t see anything on Github yet, but let me know when it’s up!

Julia has quite a few linear solving libraries, so as long as this makes matrices or LinearMaps.jl-type LinearMaps it should be fine. But what’s really needed is the interaction between meshing utilities and the discretization of operators.

16

Yeah, this sounds like the consensus as the best way forward (for now). Also, since doing this wouldn’t require too much mathematical background (or programming background, since PyCall is quite simple), this sounds like a great GSoC project. I added it to the list:

If anyone else would like to be a potential mentor for this, just submit a PR.

3 Likes
17

I recently uploaded a first version. Serious cleaning up and rationalising/documenting of the API is next on my todo list.

2 Likes
18

This is a very nice try for Julia wrapper for FEniCS, though with tiny bugs (since FEniCS is updated from time to time). I tried to debug and contacted the author.

https://github.com/JuliaDiffEq/FEniCS.jl

19

Yes, after this discussion we got a GSoC student and FEniCS.jl is what came out of it.

3 Likes
20

Has there been any progress towards a julia native FENICS like PDE solver? If not I’m finding myself in need of one, and would like to start laying down the core of an implementation. Something I’d really appreciate is some pointing in the right direction, i.e. which functionalities and Package integrations would be most desired by the community - I’d also really appreciate some pointers from people with more experience in implementing packages that take advantage of the unique opportunities offered by julia (looking at you @ChrisRackauckas) :wink:

1 Like