Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It only takes a minute to sign up.
Sign up to join this community
I've tinkered with Matlab's PDE toolbox for a while but was wondering whether anyone here has used Julia to build a PDE solver. If so, what are the advantages and limitations of Julia for PDEs? I'm mostly interested in simple PDEs with constant coefficients but am also thinking of variable coefficients for modeling purposes.
Any documentation that you could point me to would be much appreciated.
Yes, lots of people have. Automatic Jacobian sparsity handling shows up in the second tutorial of DifferentialEquations.jl, where it's able to run sparsity detection on normal Julia code to get the sparse form and perform coloring to then specialize the matrix computations. Then the tutorial ends by showing how you can swap out linear solvers for Newton-Krylov methods and add preconditioners. There's a lot more in the DifferentialEquations.jl documentation which you might find useful for time-dependent PDEs, such as GPU support, so dig around there.
Additionally there are a lot of other repos around. Oceananigans.jl is a CPU and GPU-based fluid dynamics library. CLIMA is a full climate model which uses a discrete Galerkin discretization method. NeuralNetDiffEq.jl has backwards SDE approaches for efficiently solving 1,000+ dimensional parabolic PDEs. DiffEqOperators.jl has methods for automatically constructing higher order finite difference methods with matrix-free operators. etc. Then there's things like DiffEqFlux.jl for automatically fitting models, which we used to showcase automatic discovery of PDEs.
Some of the advantages that are highlighted here is the fact that multiple dispatch and Julia's well-structured AST makes integrating compiler approaches with the solver methods quite straight forward, which means you can get arbitrary functions from the user directly in Julia code and, not only will it be fast (since it's a Julia function they're giving you), you can modify and analyze it. You can figure out its sparsity patterns, you can change their code to be more efficient for parallelism, you can reorder the equations in their code, you can replace parts by neural networks, and then throw it into traditional numerical methods. This is the kind of stuff we are doing with DifferentialEquations.jl, and I'll be putting out more resources soon on how this is being used throughout different PDE solver components and packages to improve the backend numerical methods with what are traditionally thought of as more symbolic or modeling DSL features. This is then matched with structured matrix libraries, like BlockBandedMatrices.jl, which can recompile the internals of another library to specialize the computations on not a "sparse" matrix but on structure like BlockBandedMatrices (from high dimensional PDEs), which then further improves the speed over libraries that have to treat everything as just a big unstructured sparse problem.
Those are the ones that come to mind from the top of my head. For a more comprehensive list, see the recent Survey of PDE packages.
