How composable software tools in Julia help developing multi-physics codes in geodynamics

Julia(https://julialang.org) recently emerged as a very powerful high-level computer language for (parallel) scientific computing, which allows you to “write codes like in MATLAB”, while “achieving the speed of Fortran/C”. A particular strength of Julia is that it is easy to write composable software packages that talk to each other. Here we will discuss our efforts in making Julia a development platform for geodynamic applications that significantly simplifies the process of going from a working solver to a production code which runs on massively parallel (GPU) machines.  We are working on a number of open-source packages that simplify certain steps that many geodynamics codes have in common:

• GeoParams.jl (https://github.com/JuliaGeodynamics/GeoParams.jl) is a package in which you can specify constitutive relationships (e.g., creeplaws). It automatically handles the (non-)dimensionalization of all input parameters, includes pre-defined creep laws (e.g., dislocation and diffusion creep laws), plotting routine and includes computational routines that can be directly integrated in your code.
• PETSc.jl (https://github.com/JuliaParallel/PETSc.jl) is the main interface from Julia to PETSc, including MPI support and automatic installations of PETSc (one of the main hurdles that existing users faced). We have recently extended the package to include an interface to DMSTAG, such that you create a staggered finite difference grid and assemble the stiffness matrix in a straightforward manner. You can use automatic differentiation tools in Julia to create the Jacobians for nonlinear equations, which again minimizes the required lines of code (compared to their C counterparts). At the same time, the full range of (nonlinear multigrid) PETSc solvers is available. This is thus very well suited to write implicit solvers.
• ParallelStencil.jl (https://github.com/omlins/ParallelStencil.jl) and ImplicitGlobalGrid.jl (https://github.com/eth-cscs/ImplicitGlobalGrid.jl) are packages that are devoted to solving stencils in a very efficient manner on (parallel) GPU or CPU machines, which scales to very large GPU-based computers. It is particularly efficient in combination with pseudo-transient iterative solvers and allow running codes on modern architectures.
• GeophysicalModelGenerator.jl (https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl) is a package that gives you a simple way to collect geophysical/geological data of a certain region and combine that to construct a 3D geodynamic input model setup.

Ongoing efforts include the development of a grid generation and a marker and cell advection package that work, seamlessly with both ParallelStencil and PETSc. This will allow developers to apply both direct-iterative and pseudo-transient implicit solvers to the same problem, while only having to make minimal changes to the model setup. Combined, these packages will make the step from developing a new (nonlinear) solver to having an efficient (3D) production code much easier.