CVortex.jl

A Julia wrapper for GPU accelerated vortex filament and vortex particle methods of the CVortex library.

Introduction

What does CVortex.jl do?

CVortex.jl is a wrapper for the CVortex library. It has the following functionality:

• Compute velocities induced by collections of 2D regularised vortex particles, 3D regularised vortex particles and 3D straight singular vortex filaments.
• Compute the vortex stretching term induced on 3D vortex particles by other 3D particles or vortex filaments.
• Compute the change in particle vorticity due to viscous interaction between 3D regularised vortex particles.
• Redistribute 2D and 3D vortex particles onto a grid.

What will it run on?

CVortex.jl only runs on 64bit Windows or 64bit Linux. The library is not compatible with MacOS or other CPU instruction sets.

To obtain the maximum benefit you'll also need an OpenCL 1.2 compatible GPU or iGPU. This includes:

• Intel integrated graphics
• AMD integrated graphics and discrete GPUs
• Nvidia GPUs (Any GPU that runs CUDA)

You'll also have needed to have installed the appropriate GPU drivers. Note that many hypervisors (programs that allow you to run virtual machines such as VirtualBox) don't pass through graphics hardware.

Even if you don't have compatable a compatible GPU, you'll still benefit from the multicore implementation.

How can I get CVortex.jl?

You'll need to add the package to your system. Run Julia and:

(v1.1) pkg> add CVortex

remembering that to access to package environment the ] must be used. Binaries for the CVortex library will automatically be downloaded.

Is there documentation?

Yes! The ordinary help syntax within Julia (Type ? within the REPL) will give you help on a particular topic. For example:

help?> particle_induced_velocity

For examples see CVortex.jl examples

Using CVortex.jl

The first thing you'll need to do is to import CVortex.jl. This can be done using

using CVortex

Vortex filaments

All vortex filaments in CVortex are straight and singular. They have three properties, a start point, an end point and a vorticity. The first two are 3D, and the latter is a scalar.

Velocity

To obtain the velocity induced on a point in 3D one uses:

startp = rand(3)  # Filament start coordinate
endp = rand(3)    # Filament end coordinate
fvort = rand()    # A scalar. Filament's vorticity
mesp = rand(3)    # Velocity measurement location.
vel = filament_induced_velocity(startp, endp, fvort, mesp);

The returned vel is a Float32 array of length 3.

We use our hardware better if we group computations together. Suppose we had N vortex filaments, we can vectorise the computation of the influence on a measurement point as

N = 10000
startps = rand(N, 3)
endps = rand(N, 3)
fvorts = rand(N)
mesp = rand(3)
vel = filament_induced_velocity(startps, endps, fvorts, mesp)

Again, vel is a Float32 array of length 3.

To create a problem suitable for GPU accelleration, we need multiple-multiple problems. To do this the measurement points becomes a matrix:

N = 10000
M = 100000
startps = rand(N, 3)
endps = rand(N, 3)
fvorts = rand(N)
mesp = rand(M, 3)
vels = filament_induced_velocity(startps, endps, fvorts, mesp)

where vels is an M by 3 Float32 matrix.

Influence matrix

Its often desirable to obtain the influence of a vortex filament (perhaps it the context of a vortex ring) on the velocity in a given direction at a point. A function for this is included:

nvels = induced_velocity_influence_matrix(
      filament_start_coords :: Matrix{<:Real},
      filament_end_coords :: Matrix{<:Real},
      measurement_points :: Matrix{<:Real},
      measurement_directions :: Matrix{<:Real})

For N vortex filaments, a matrix A is returned such that, for a length N vector of filament vorticity called b, the induced velocities measured at the M points in M directions would be given by A*b.

Vortex particles

Vortex particles are blobs of vorticity. Whilst they can be singular, this isn't good for long term stability. CVortex therefore implements vortex particle regularisation.

Regularisation

CVortex.jl uses the struct RegularisationFunction to group together functions relevent to a regularisation method. A RegularisationFunction can be obtained using pre-programmed routines:

singular_reg = singular_regularisation()
planet_reg = planetary_regularisation()
winckel_reg = winckelmans_regularisation()
gauss_reg = gaussian_regularisation()

singular_regularisation() isn't actually a regularisation because it allows the use of singular vortex particles. planetary_regularisation() allow regularisation, but cannot be used in viscous schemes. winckelmans_regularisation() is a higher order algebraic regularisation. gaussian_regularisation() is normal gaussian regularisation.

2D vortex particles

Using 2D vortex particles is much like using singular vortex filament, with two additions:

• A regularisation function is required
• A regularisation distance is required

The regularisation distance is like the radius of the vortex particles. Roughly, it represents the finest fidelity that the field can resolve. Consequently, vortex particles must overlap for a good solution.

Velocity

To obtain the velocity induced by a 2D vortex particle one uses

vels = particle_induced_velocity(particle_positions, particle_vorts, 
  measurement_points, regularisation, regularisation_distance)

Where:

• For single particle -> single measurement point: particle_positions is a length 2 vector, particle_vorts is scalar and measurement_points is a length 2 vector. vels is a length 2 vector.
• For multiple particles -> single measurement point: particle_positions is an N by 2 matrix, particle_vorts is a length N vector and measurement_points is a length 2 vector. vels is a length 2 vector.
• For multiple particles -> multiple measurement points: particle_positions is an N by 2 matrix, particle_vorts is a length N vector and measurement_points is an M by 2 matrix. vels is an M by 2 matrix. In all cases, regularisation_distance is as scalar. Different sized particles aren't supported.

Viscous rate of change of vorticity

For viscous vortex particle method, the rate of change of vorticity of each particle is computed.

dvorts = particle_visc_induced_dvort(
  inducing_particle_positions, inducing_particle_vorts, inducing_particle_areas,
  induced_particle_positions, induced_particle_vorts, induced_particle_areas,
  regularisation, regularisation_distance, kinematic_viscosity)

Here the rate of change of vorticity on the "induced" particles is computed. The particle_area variables is of the same type as the vorticity vector.

3D Vortex particles

3D vortex particles are characterised by a position (in 3D) and a vorticity vector (again, in 3D). Like 2D particles, a regularisation function and a regularisation distance is required for computation.

Velocity

To obtain the velocity induced by a 3D vortex particle one uses

vels = particle_induced_velocity(particle_positions, particle_vorts, 
  measurement_points, regularisation, regularisation_distance)

Where:

• For single particle -> single measurement point: particle_positions is a length 3 vector, particle_vorts is a length 3 vector and measurement_points is a length 3 vector. vels is a length 3 vector.
• For multiple particles -> single measurement point: particle_positions is an N by 3 matrix, particle_vorts is an N by 3 matrix and measurement_points is a length 3 vector. vels is a length 3 vector.
• For multiple particles -> multiple measurement points: particle_positions is an N by 3 matrix, particle_vorts is an N by 3 matrix and measurement_points is an M by 3 matrix. vels is an M by 3 matrix. In all cases, regularisation_distance is as scalar. Different sized particles aren't supported.

Vortex stretching

In 3D vorticies can be "stretched" leading to a rate of change of vorticity. To compute this the following is used:

dvort = particle_induced_dvort(
        inducing_particle_position, inducing_particle_vorticity,
        induced_particle_position, induced_particle_vorticity,
        kernel :: RegularisationFunction, regularisation_radius :: Real)

Viscous rate of change of vorticity

For viscous vortex particle method, the rate of change of vorticity of each particle is computed.

dvorts = particle_visc_induced_dvort(
  inducing_particle_positions, inducing_particle_vorts, inducing_particle_areas,
  induced_particle_positions, induced_particle_vorts, induced_particle_areas,
  regularisation, regularisation_distance, kinematic_viscosity)

Here the rate of change of vorticity on the "induced" particles is computed. The particle_area variables is of the same type as the vorticity vector.

Vortex particle redistribution

Vortex particles can be redistributed onto a grid to fix problems introduced by particles spreading out of grouping together.

To do this some kind of redistribution scheme is required. This scheme is encapsulated within a RedistributionFunction struct. These can be created as

scheme = lambda0_redistribution()
scheme = lambda1_redistribution()
scheme = lambda2_redistribution()
scheme = lambda3_redistribution()
scheme = m4p_redistribution()

The lambda3_redistribution() scheme and m4p_redistribution() scheme are generally recommended. The lambda0_redistribution() and lambda2_redistribution() are discontinious, and so can cause problems numerically. The lambda1_redistribution() is dissipative.

Having chosen a scheme, particles can then be redistributed:

positions, vorts, areas = redistribute_particles_on_grid(
        particle_positions, inducing_particle_vorticity,
        redistribution_function, grid_density;
        negligible_vort=1e-4, max_new_particles=-1)

There are two optional parameters, negligible_vort and max_new_particles. These are designed to stop lots of vortex particles with very small vorticities being created.

negligible_vort is a threshold for discarding particles, and should be a value between 0 (discard no particles) and 1 (discard all particles). It is implemented as the proportion of the average particle's vorticity that any particle must have possess to be kept. The vorticity of any discarded particle is distributed evenly among the remaining particles.

max_new_particles is a hard limit on the number of new vortex particles that can be created. When equal to -1, there is no limit. If negligible_vort is chosen such that there are more than the max_new_particles remaining, further particles are discarded until the number of particles is less than max_new_particles. Due to the implementation, this may result in fewer particles than max_new_particles.

Mixing 3D vortex particles and filaments

It is possible to mix vortex particles and filaments in 3D. Since vortex filaments are singular, it is not possible to include viscosity for these problems (unless you're willing to cheat somehow).

The only addition function required for putting both in one problem is the vortex stretching induced by vortex filaments on vortex particles. This can be obtained using

dvort = filament_induced_dvort(
      filament_start_coord,
      filament_end_coord,
      filament_strength,
      induced_particle_position ,
      induced_particle_vorticity)

where everything is assumed to be singular.

Contolling the accelerators / GPUs

In computers with multiple GPUs (probably an iGPU + discrete GPU) it may be desirable to control which GPU is being used, or in some cases to stop GPUs being used at all.

But first, one must know how many GPUs CVortex has found:

number_of_gpus = number_of_accelerators()

where an integer is returned.

The accelerators are given an index of 1:number_of_accelerators(). Acclerators can then be controlled and investigated using the index.

To obtain the name one uses:

name = accelerator_name(accelerator_index)

Note that the name may not be unique among your GPUs, or even share the name of the product you purchased. For example, an for an AMD RX Vega 56:

julia> accelerator_name(1)
"gfx900"

To investigate whether CVortex is using a GPU:

in_use = accelerator_enabled(index)

which returns 1 (true) or 0 (false).

To enable an accelerator

accelerator_enable(index)

and to disable an accelerator

accelerator_disable(index)

Potentially FAQ

Why does CVortex.jl return Float32s? Because (almost) all the underlying code uses floats. GPU manufactures cripple the double precision speed of their consumer-level GPUs, or may not include double-precision capability at all (it isn't required by the OpenCL spec.). But since the discretisation of vortex particles of vortex filaments lead to more error than single precision computing, the cost of using single precision is negligible, whilst the compatability/performance benifits are huge.

Why is the program hanging? The implementation of some OpenCL drivers is supposedly dodgy, especially on older GPUs.

Why isn't CVortex.jl available on MacOS or X architecture? CVortex is theoretically compatable with MacOS and any CPU architecture that can be targetet by MSVC, GCC or Clang. For multithreading, OpenMP has to be available. For GPU acceleration, an implementation of the full OpenCL 1.2 profile has to be available. CVortex isn't compiled for these architectures or Mac, but I lack the hardware on which to compile or test binaries. For a sufficiently keen bean, it would be possible to compile cvortex from source, and copy the shared library into CVortex.jl build directory.

Why isn't it faster? It would be possible to make CVortex faster, but to do so would complicate the interface. The code has also not be tailored to any particular hardware.

Why does using CVortex take such a long time to run? using CVortex calls the underlying CVortex library's initialisation function which must in turn compile the OpenCL kernels for used by the programs.

Why isn't CVortex using my iGPU/GPU Thats potentially a big questions. If number_of_accelerators() returns the expected number of devices, its probably because the problem isn't suitable for GPU acceleration (only multiple-multiple problems are accelerated, and even then the problem must have a sufficient number of both target and source objects). If the number is less than expected, things get more complicated. It may be that the OpenCL kernels could not be sucessfully compiled (I've not encountered this) or, more likely, the OpenCL ICD loader didn't find the device's OpenCL runtime library. This might be because the drivers aren't properly installed. Additionally, on Windows, driver installers are liable to overwrite files installed by other driver installers. If you're running a virtual machine, check that the GPU is being passed through (if the hypervisor is even capable of doing this).