# Achieving C++ Speeds for a CFD Algorithm: Julia spends too much time on pre-allocation?

**URL:** <https://discourse.julialang.org/t/achieving-c-speeds-for-a-cfd-algorithm-julia-spends-too-much-time-on-pre-allocation/17410>\
**Category:** Performance\
**Created:** [November 11, 2018, 7:51pm UTC](https://discourse.julialang.org/t/achieving-c-speeds-for-a-cfd-algorithm-julia-spends-too-much-time-on-pre-allocation/17410 "2018-11-11T19:51:55Z")\
**Posts on this page:** 3\
**Page:** 2

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [November 12, 2018, 4:04am UTC](https://discourse.julialang.org/t/achieving-c-speeds-for-a-cfd-algorithm-julia-spends-too-much-time-on-pre-allocation/17410/21 "2018-11-12T04:04:51Z")

</div>

> [@rdeits](#):
>
> For example, running Julia with `julia -O3` can enable additional SIMD optimizations ([source](https://docs.julialang.org/en/v1.0/base/simd-types/#SIMD-Support-1))

That documentation is out of date and the comment about `O3` removed. The SLP optimizer pass is enabled by default now.

---

<div class="post-metadata">

**Author:** ![EmadMasroor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/emadmasroor/32/7527_2.png) [@EmadMasroor](https://discourse.julialang.org/u/EmadMasroor)\
**Post date:** [November 28, 2018, 2:33am UTC](https://discourse.julialang.org/t/achieving-c-speeds-for-a-cfd-algorithm-julia-spends-too-much-time-on-pre-allocation/17410/22 "2018-11-28T02:33:38Z")

</div>

Hi -

Thanks to those of you who helped me make sense of the compilation/pre-allocation questions I had earlier. Once I found some time to make changes to the flow of my program, I have come up with the following logical flow.

The following code is in a file called `main.jl`

```julia
include("GaussSiedel.jl"); #
include("functions.jl"); # 
include("Loop.jl"); # these are files which contain all of my function declarations

# The following are values that don't change.
const delx = 1/64;
const dely = 1/64;
const x = 0:delx:1;
const y = 0:dely:1;
const dV = delx*dely;
const Re = 100;
const Co = 0.9;
const ncol = length(x);
const nrow = length(y);
const nc = nrow*ncol;
const interiorCells = nc - (2*nrow + 2*ncol - 4);
# In the following line, I initialize and pre-allocate all of the arrays in the problem.
u,v,oldu,oldv,deltau,deltav,Convx,Convy,Diffx,Diffy,Divergence,p,residue,gradxP,gradyP = initializeArrays(delx,dely);

# Here, I declare the "main" function for my program.
function main(u,v,oldu,oldv,deltau,deltav,Convx,Convy,Diffx,Diffy,Divergence,p,residue,gradxP,gradyP)
deltaVL2Norm = 1;
EnforceBoundaryConditions(u,v,p,nrow,ncol,nc);
coefficients = stencilFunction(delx,dely);
# Here is the while loop in which my function "Loop" is called repeatedly. All substantative calculations happen in the Loop function.
while deltaVL2Norm > 1e-8
	deltaVL2Norm = Loop(deltaVL2Norm,u,v,oldu,oldv,deltau,deltav,Convx,Convy,Diffx,Diffy,Divergence,p,residue,gradxP,gradyP);
end
printResults(u,v,p,nrow,ncol);
end

```

The way I run my code is by first typing `include("main.jl")` at the Julia REPL. This, I believe, pre-allocates arrays, defines constants, and (crucially) declares the function `main(...)` which in turn declares `Loop(...)`. Next, I type `main()` together with all of its arguments at the Julia REPL, which executes everything.

Currently, this method gives me fairly good performance and accurate results compared to a benchmark CFD problem for coarse grid sizes. For finer grid sizes, the iterative solver takes too many iterations. While this is most probably my own fault, it is certainly taking much, much more time than C++ does for this problem.

I have the following questions:

- When I execute `main(...)` once, the function is being executed for the first time and so it is also compiled at this time. However, I can’t execute `main(...)` again because now, all my variables are already at their correct values and the `Loop()` function will only be called once, since all variables satisfy the physical solution to this problem. One solution to this could be to place the initializing line (which resets all variables to zero)  
`u,v,oldu,oldv,deltau,deltav,Convx,Convy,Diffx,Diffy,Divergence,p,residue,gradxP,gradyP = initializeArrays(delx,dely);`  
inside the `main()` function. This way, when I call `main(...)` a second time, it will be re-solving the problem from scratch, and it should not be being compiled and so should be **much** faster. BUT when I do this, the following things happen. Why?

1. It takes the same time in the second try as in the first, so my intuition about compiling of `main(...)` is wrong.
2. Curiously, the solution is weird. The Poisson solver for pressure (an iterative Gauss-Siedel method that occurs in each time step) no longer converges in 2 iterations toward the end of the algorithm (when things are reaching steady state and the pressure equation should be near-immediately satisfied, physically) and instead always takes at least 28 iterations. This is extremely weird and I cannot see why placing the function `initializeArrays(delx,dely)` inside `main(...)` should make any difference.

- It appears to make no difference to my code whether or not the arrays `u,v,oldu,oldv` etc. are included as arguments to `main(..)` or not. Even when I define `main()` to have no arguments, the code works exactly the same. What’s going on? Are my arrays (which are the crucial variables which the `Loop()` function operates on and modifies at each time step) living in the ‘right’ place, or are they floating around in relatively inaccessible parts of the memory?

- Does it make sense to define constants using the keyword `const` in my program?

- I am finding it harder and harder to stick to my original goal with this project, which (as the title of the original post states) was to achieve C++ speeds for an identical algorithm in Julia, at least once the code has been compiled. However, as written, there appears to be no good way to do this if I put everything into functions, as I believe I should.

Any other advice about how to improve my code, make it more ‘Julian’, and improve things like memory management would be highly appreciated.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 28, 2018, 6:14am UTC](https://discourse.julialang.org/t/achieving-c-speeds-for-a-cfd-algorithm-julia-spends-too-much-time-on-pre-allocation/17410/23 "2018-11-28T06:14:18Z")

</div>

- It does make sense to define your `const`s
- Your `deltaVL2Norm` should probably start out as a `Float64`, i.e., `deltaVL2Norm = 1.0`
- If `Loop` takes a long time to run, the dynamic dispatch that happens due to lookup of the global variables in the arguments is probably insignificant in comparison to the runtime of `Loop`, hence you might not see a speedup when sending the arrays as arguments to `main`.
- I can not see why it should make any difference initializing the arrays inside `main`, I have no idea what the function `initializeArrays` does.

[Previous page](https://discourse.julialang.org/t/achieving-c-speeds-for-a-cfd-algorithm-julia-spends-too-much-time-on-pre-allocation/17410.md?page=1)