Adventures in Julia Programming

So some early testing of a better optimised GasCode program suggests that code like: ∫▽●V[i] = 0.0 Is really badly optimised. Instead, you're better off using ∫▽●V = zeros(ncells) Or potentially @inbounds @sync @distributed for i in 1:ncells ∫▽●V[i] = 0.0 end If you're using the Distributed package Note that if you don't include the @sync in this style of distributed loop, you don't actually update the array.

Let's scale things back, and look at a simple example - brute force calculation of Pi. We can do this using the following formula: This is easy enough to code up in various languages: C++ #include <iostream> #include <stdio.h> static long nsteps = 10000000000; int main() { long i; double x, pisum = 0.0; double mypi; double pi = 3.141592653589793; double step = 1.0/(double)nsteps; for (i=0;i<nsteps;i++) { x = (i+0.5)*step; pisum = pisum + 4.0/(1.0+x*x); } mypi = pisum * step; printf("%20.16f\n", mypi); printf("%20.16e\n", mypi - pi); } Fortran program calculate_pi use, intrinsic :: iso_fortran_env USE OMP_LIB implicit none integer(kind=INT64), parameter :: nsteps = 10000000000_INT64 integer(kind=INT64) :: i real(kind=REAL64) :: stepsize, x, pisum real(kind=REAL64) :: t1, t2 real(kind=REAL64) :: mypi, diff real(kind=REAL64) :: pi = 3.141...

So, we have basic, Lagrangian-only versions of GasCode written in Fortran, Python and Julia, so we can take a look at their relative performance. For these comparisons, I'll use the spherical Sod problem on a 50x50 mesh: This problem is being used simply because it's the hard-coded one, and I've not added a mesh reader to the Python or Fortran versions of the code yet. Input options are as follows: # EoS gamma = 1.4 # Artificial Viscosity CL = 0.5 CQ = 0.75 # Mesh Extents x0 = 0.0 x1 = 1.0 y0 = 0.0 y1 = 1.0 # Timestep, etc. t0 = 0.0 tend = 0.205 dtinit = 0.0001 dtmax = 0.01 growth = 1.02 # Mesh resolution nregions = 1 meshx = 50 meshy = 50 # Debug settings debug_step_count = 0   # Set to zero to disable debugging # Cutoffs rhocutoff = 1.0e-6 I.e. basic monotonic Q, and no anti-hourglass. Using these settings, and running the Fortran code along with the most simple (/naive) versions of the Julia and Python codes, we get the following runtimes for just the Lagrangian ...

So for ALE/SALE, we need a connectivity array linking nodes together. Let's look again at part of our mesh: This time, I have labelled both the cells and the nodes. Let's consider node 114 - we need to create a node-node connectivity, nodnodconn, such that: nodnodconn[1, 114] = 113 nodnodconn[2, 114] = 118 nodnodconn[3, 114] = 56 nodnodconn[4, 114] = 57 We currently have two pieces of information which we can use to construct this connectivity: The identities of the nodes surrounding any given cell (nodelist) The orthogonal neighbours of any given cell (elelconn)  Let's look at how we might map our nodnodconn array onto a simple 4 cell/9 node stencil: Here, we would be looping over cells rather than nodes, and our aim is to create the nodnodconn entries for N0. I.e: nodnodconn[1, N0] = N1 nodnodconn[2, N0] = N2 nodnodconn[3, N0 ] = N3 nodnodconn[4, N0 ] = N4 The positions of the node labels in the diagram shows the easiest way to access them given ...

I did this yesterday, but only started the blog today... I t is possible to create an executable file from a Julia script, with a little bit of effort. Basically, we need to create a Julia package project, and then use the   PackageCompiler   library to turn this into an executable:   Process to create a Julia package Launch an interactive julia session Press the "]" key, the Julia prompt will be replaced with a pkg prompt: ( @1.4 ) pkg> Type "generate <packagename>" This will create a new package with the given name. In doing so, it will create a directory with the same name, which will include a "src" subdirectory, and a Project.toml file which describes the project to Julia Activate the project by typing "activate <packagename>" Add any and all dependencies used by the project. This is done via the "add <package>" command, which can accept more than one package name at a time. This will update ...

So for certain features of the hydrodynamic solve, such as a Christensen Q, and ALE/SALE, we need some form of element-element connectivity array. Consider the following mesh fragment: The aim of the connectivity array is to link cell 13 to its four orthogonal neighbours, 8, 12, 14 and 18. Specifically, however, we want an array of size [4,ncells] where the first index corresponds to the directions down, right, up, left in that order. I.e: elelconn[1,13] = 12 elelconn[2,13] = 18 elelconn[3,13] = 14 elelconn[4,13] = 8 In order to calculate this array correctly, we first need to look at how GMSH generates the cell-node connectivity. Turns out we have control of this, via the Transfinite Surface input. We'll set up the input file so that it creates the following node ordering (in nodelist): Here the green numbers show how the nodes are ordered (i.e. indices of nodelist), whilst the red numbers show (again) how we want the edges to be numbered (i.e. indices of elelcon...