More examples
The examples folder contains further examples, each illustrating a particular feature. A simple boolean USE_GPU defines whether they run on GPU or CPU. Most of them report the effective memory throughput T_eff (see Performance metric), and the multi-xPU examples rely on ImplicitGlobalGrid.jl.
Memory copy
Memory copy benchmark
A memory copy kernel written with @parallel and math-close notation, measuring the effective memory throughput T_eff reachable on the hardware (a practical upper bound for stencil computations): memcopy3D.jl.
Code
julia
const USE_GPU = true
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
@init_parallel_stencil(CUDA, Float64, 3);
else
@init_parallel_stencil(Threads, Float64, 3);
end
@parallel function copy3D!(T2, T, Ci)
@all(T2) = @all(T) + @all(Ci);
return
end
function memcopy3D()
# Numerics
nx, ny, nz = 512, 512, 512; # Number of gridpoints in dimensions x, y and z
nt = 100; # Number of time steps
# Array initializations
T = @zeros(nx, ny, nz);
T2 = @zeros(nx, ny, nz);
Ci = @zeros(nx, ny, nz);
# Initial conditions
Ci .= 1/2.0;
T .= 1.7;
T2 .= T;
# Time loop
for it = 1:nt
if (it == 11) global t0=time(); end # Start measuring time.
@parallel copy3D!(T2, T, Ci);
T, T2 = T2, T;
end
time_s=time()-t0
# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(Data.Number); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = time_s/(nt-10); # Execution time per iteration [s]
T_eff = A_eff/t_it; # Effective memory throughput [GB/s]
println("time_s=$time_s T_eff=$T_eff");
end
memcopy3D()Memory copy benchmark with @parallel_indices
The same benchmark written with @parallel_indices: memcopy3D_parindices.jl.
Code
julia
const USE_GPU = true
using ParallelStencil
@static if USE_GPU
@init_parallel_stencil(CUDA, Float64, 3);
else
@init_parallel_stencil(Threads, Float64, 3);
end
@parallel_indices (ix,iy,iz) function copy3D!(T2, T, Ci)
T2[ix,iy,iz] = T[ix,iy,iz] + Ci[ix,iy,iz];
return
end
function memcopy3D()
# Numerics
nx, ny, nz = 512, 512, 512; # Number of gridpoints in dimensions x, y and z
nt = 100; # Number of time steps
# Array initializations
T = @zeros(nx, ny, nz);
T2 = @zeros(nx, ny, nz);
Ci = @zeros(nx, ny, nz);
# Initial conditions
Ci .= 1/2.0;
T .= 1.7;
T2 .= T;
# Time loop
for it = 1:nt
if (it == 11) global t0=time(); end # Start measuring time.
@parallel copy3D!(T2, T, Ci);
T, T2 = T2, T;
end
time_s=time()-t0
# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(Data.Number); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = time_s/(nt-10); # Execution time per iteration [s]
T_eff = A_eff/t_it; # Effective memory throughput [GB/s]
println("time_s=$time_s T_eff=$T_eff");
end
memcopy3D()Logical arrays of small arrays / structs
A memory copy of CellArrays allocated with the keyword celldims and cell types defined with @CellType, using math-close notation as well as explicit indexing of cells and of cell fields: memcopyCellArray3D.jl.
Code
julia
const USE_GPU = true
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
@init_parallel_stencil(CUDA, Float64, 3);
else
@init_parallel_stencil(Threads, Float64, 3);
end
@parallel function copy3D!(T2::Data.CellArray, T::Data.CellArray, Ci::Data.CellArray)
@all(T2) = @all(T) + @all(Ci);
return
end
@parallel_indices (ix,iy,iz) function copy3D_explicit!(T2::Data.CellArray, T::Data.CellArray, Ci::Data.CellArray)
T2[ix,iy,iz] = T[ix,iy,iz] + Ci[ix,iy,iz]
return
end
using CellArrays
@parallel_indices (ix,iy,iz) function copy3D_large_cells!(T2::Data.CellArray, T::Data.CellArray, Ci::Data.CellArray)
for cy = 1:cellsize(T2,2), cx = 1:cellsize(T2,1)
field(T2,cx,cy)[ix,iy,iz] = field(T,cx,cy)[ix,iy,iz] + field(Ci,cx,cy)[ix,iy,iz]
end
return
end
@CellType Cell4x4 fieldnames=(xxxx, yxxx, xyxx, yyxx, xxyx, yxyx, xyyx, yyyx, xxxy, yxxy, xyxy, yyxy, xxyy, yxyy, xyyy, yyyy) dims=(4,4)
@CellType TCell4x4 fieldnames=(xxxx, yxxx, xyxx, yyxx, xxyx, yxyx, xyyx, yyyx, xxxy, yxxy, xyxy, yyxy, xxyy, yxyy, xyyy, yyyy) dims=(4,4) parametric=true
function memcopy3D()
# Numerics
nx, ny, nz = 128,128,1024 #CPU:128,128,128 # Number of gridpoints in dimensions x, y and z
celldims = (4,4)
nt = 100; # Number of time steps
# Array initializations
T = @zeros(nx, ny, nz, celltype=Cell4x4); # or: celltype=TCell4x4{Float64}
T2 = @zeros(nx, ny, nz, celltype=Cell4x4); # or: celltype=TCell4x4{Float64}
Ci = @zeros(nx, ny, nz, celldims=celldims);
# Initial conditions
@fill!(Ci, 0.5);
@fill!(T, 1.7);
copy!(T2.data, T.data);
# Time loop
for it = 1:nt
if (it == 11) global t0=time(); end # Start measuring time.
@parallel copy3D!(T2, T, Ci); # or: @parallel copy3D_explicit!(T2, T, Ci); or: @parallel copy3D_large_cells!(T2, T, Ci);
T, T2 = T2, T;
end
time_s=time()-t0
# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*nz*prod(celldims)*sizeof(Data.Number); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = time_s/(nt-10); # Execution time per iteration [s]
T_eff = A_eff/t_it; # Effective memory throughput [GB/s]
println("time_s=$time_s T_eff=$T_eff");
end
memcopy3D()Named tuples of arrays
A kernel operating on named tuples of arrays (Data.NamedArrayTuple), e.g. on the components of a vector field: multimemcopyNamedArrayTuple3D.jl.
Code
julia
const USE_GPU = true
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
@init_parallel_stencil(CUDA, Float64, 3);
else
@init_parallel_stencil(Threads, Float64, 3);
end
@parallel function copy3D!(V2::Data.NamedArrayTuple, V::Data.NamedArrayTuple, dV::Data.NamedArrayTuple)
@all(V2.x) = @all(V.x) + @all(dV.x)
@all(V2.y) = @all(V.y) + @all(dV.y)
@all(V2.z) = @all(V.z) + @all(dV.z)
return
end
function memcopy3D()
# Numerics
nx, ny, nz = 512, 512, 512; # Number of gridpoints in dimensions x, y and z
nt = 100; # Number of time steps
# Array initializations
V = (x = @zeros(nx+1, ny, nz), y = @zeros(nx, ny+1, nz), z = @zeros(nx, ny, nz+1))
V2 = (x = @zeros(nx+1, ny, nz), y = @zeros(nx, ny+1, nz), z = @zeros(nx, ny, nz+1))
dV = (x = @zeros(nx+1, ny, nz), y = @zeros(nx, ny+1, nz), z = @zeros(nx, ny, nz+1))
# Initial conditions
dV.x .= 1/2.0; dV.y .= 1/2.0; dV.z .= 1/2.0;
V.x .= 1.7; V.y .= 1.7; V.z .= 1.7;
V2.x .= V.x; V2.y .= V.y; V2.z .= V.z;
# Time loop
for it = 1:nt
if (it == 11) global t0=time(); end # Start measuring time.
@parallel copy3D!(V2, V, dV);
V, V2 = V2, V;
end
time_s=time()-t0
# Performance
A_eff = (2*3+3)*1/1e9*nx*ny*nz*sizeof(Data.Number); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = time_s/(nt-10); # Execution time per iteration [s]
T_eff = A_eff/t_it; # Effective memory throughput [GB/s]
println("time_s=$time_s T_eff=$T_eff");
end
memcopy3D()Heat diffusion
Shared memory
A 2-D heat diffusion solver with a @parallel_indices kernel that explicitly uses shared memory (@sharedMem): diffusion2D_shmem_novis.jl.
Code
julia
const USE_GPU = true
using BenchmarkTools
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
@init_parallel_stencil(CUDA, Float64, 2);
else
@init_parallel_stencil(Threads, Float64, 2);
end
@parallel_indices (ix,iy) function diffusion2D_step!(T2, T, Ci, lam, dt, _dx, _dy)
tx = @threadIdx().x + 1
ty = @threadIdx().y + 1
T_l = @sharedMem(eltype(T), (@blockDim().x+2, @blockDim().y+2))
T_l[tx,ty] = T[ix,iy]
if (1<ix<size(T2,1) && 1<iy<size(T2,2))
if (@threadIdx().x == 1) T_l[tx-1,ty] = T[ix-1,iy] end
if (@threadIdx().x == @blockDim().x) T_l[tx+1,ty] = T[ix+1,iy] end
if (@threadIdx().y == 1) T_l[tx,ty-1] = T[ix,iy-1] end
if (@threadIdx().y == @blockDim().y) T_l[tx,ty+1] = T[ix,iy+1] end
@sync_threads()
T2[ix,iy] = T_l[tx,ty] + dt*(Ci[ix,iy]*(
- ((-lam*(T_l[tx+1,ty] - T_l[tx,ty])*_dx) - (-lam*(T_l[tx,ty] - T_l[tx-1,ty])*_dx))*_dx
- ((-lam*(T_l[tx,ty+1] - T_l[tx,ty])*_dy) - (-lam*(T_l[tx,ty] - T_l[tx,ty-1])*_dy))*_dy
));
end
return
end
function diffusion2D()
# Physics
lam = 1.0; # Thermal conductivity
c0 = 2.0; # Heat capacity
lx, ly = 1.0, 1.0; # Length of computational domain in dimension x and y
# Numerics
nx, ny = 512*32, 512*32; # Number of gridpoints in dimensions x and y
nt = 100; # Number of time steps
dx = lx/(nx-1); # Space step in x-dimension
dy = ly/(ny-1); # Space step in y-dimension
_dx, _dy = 1.0/dx, 1.0/dy;
# Array initializations
T = @zeros(nx, ny);
T2 = @zeros(nx, ny);
Ci = @zeros(nx, ny);
# Initial conditions
Ci .= 1/c0; # 1/Heat capacity
T .= 1.7;
T2 .= T; # Assign also T2 to get correct boundary conditions.
#
# GPU launch parameters
threads = (32, 8)
blocks = (nx, ny) .÷ threads
#
# Time loop
dt = min(dx^2,dy^2)/lam/maximum(Ci)/4.1; # Time step for 2D Heat diffusion
for it = 1:nt
if (it == 11) GC.enable(false); global t_tic=time(); end # Start measuring time.
@parallel blocks threads shmem=prod(threads.+2)*sizeof(Float64) diffusion2D_step!(T2, T, Ci, lam, dt, _dx, _dy);
T, T2 = T2, T;
end
time_s = time() - t_tic
# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*sizeof(eltype(T)); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = time_s/(nt-10); # Execution time per iteration [s]
T_eff = A_eff/t_it; # Effective memory throughput [GB/s]
println("time_s=$time_s t_it=$t_it T_eff=$T_eff");
# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*sizeof(eltype(T)); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = @belapsed begin @parallel $blocks $threads shmem=prod($threads.+2)*sizeof(Float64) diffusion2D_step!($T2, $T, $Ci, $lam, $dt, $_dx, $_dy); end
println("Benchmarktools (min): t_it=$t_it T_eff=$(A_eff/t_it)");
end
diffusion2D()Multi-xPU with performance measurement
The 50-lines multi-xPU example with a measurement of T_eff: diffusion3D_multigpucpu_novis.jl.
Code
julia
const USE_GPU = true
using ImplicitGlobalGrid
import MPI
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
@init_parallel_stencil(CUDA, Float64, 3);
else
@init_parallel_stencil(Threads, Float64, 3);
end
@parallel function diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz)
@inn(T2) = @inn(T) + dt*(lam*@inn(Ci)*(@d2_xi(T)*_dx^2 + @d2_yi(T)*_dy^2 + @d2_zi(T)*_dz^2));
return
end
function diffusion3D()
# Physics
lam = 1.0; # Thermal conductivity
c0 = 2.0; # Heat capacity
lx, ly, lz = 1.0, 1.0, 1.0; # Length of computational domain in dimension x, y and z
# Numerics
nx, ny, nz = 512, 512, 512; # Number of gridpoints in dimensions x, y and z
nt = 100; # Number of time steps
me, dims = init_global_grid(nx, ny, nz);
dx = lx/(nx_g()-1); # Space step in x-dimension
dy = ly/(ny_g()-1); # Space step in y-dimension
dz = lz/(nz_g()-1); # Space step in z-dimension
_dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz;
# Array initializations
T = @zeros(nx, ny, nz);
T2 = @zeros(nx, ny, nz);
Ci = @zeros(nx, ny, nz);
# Initial conditions
Ci .= 1/c0; # 1/Heat capacity
T .= 1.7;
T2 .= T; # Assign also T2 to get correct boundary conditions.
# Time loop
dt = min(dx^2,dy^2,dz^2)/lam/maximum(Ci)/8.1; # Time step for 3D Heat diffusion
for it = 1:nt
if (it == 11) tic(); end # Start measuring time.
@parallel diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz);
update_halo!(T2);
T, T2 = T2, T;
end
time_s = toc()
# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(Data.Number); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = time_s/(nt-10); # Execution time per iteration [s]
T_eff = A_eff/t_it; # Effective memory throughput [GB/s]
if (me==0) println("time_s=$time_s T_eff=$T_eff"); end
finalize_global_grid();
end
diffusion3D()Hiding communication
The communication is hidden behind computation with @hide_communication: diffusion3D_multigpucpu_hidecomm_novis.jl.
Code
julia
const USE_GPU = true
using ImplicitGlobalGrid
import MPI
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
@init_parallel_stencil(CUDA, Float64, 3);
else
@init_parallel_stencil(Threads, Float64, 3);
end
@parallel function diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz)
@inn(T2) = @inn(T) + dt*(lam*@inn(Ci)*(@d2_xi(T)*_dx^2 + @d2_yi(T)*_dy^2 + @d2_zi(T)*_dz^2));
return
end
function diffusion3D()
# Physics
lam = 1.0; # Thermal conductivity
c0 = 2.0; # Heat capacity
lx, ly, lz = 1.0, 1.0, 1.0; # Length of computational domain in dimension x, y and z
# Numerics
nx, ny, nz = 512, 512, 512; # Number of gridpoints in dimensions x, y and z
nt = 100; # Number of time steps
me, dims = init_global_grid(nx, ny, nz);
dx = lx/(nx_g()-1); # Space step in x-dimension
dy = ly/(ny_g()-1); # Space step in y-dimension
dz = lz/(nz_g()-1); # Space step in z-dimension
_dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz;
# Array initializations
T = @zeros(nx, ny, nz);
T2 = @zeros(nx, ny, nz);
Ci = @zeros(nx, ny, nz);
# Initial conditions
Ci .= 1/c0; # 1/Heat capacity
T .= 1.7;
T2 .= T; # Assign also T2 to get correct boundary conditions.
# Time loop
dt = min(dx^2,dy^2,dz^2)/lam/maximum(Ci)/8.1; # Time step for 3D Heat diffusion
for it = 1:nt
if (it == 11) tic(); end # Start measuring time.
@hide_communication (16, 2, 2) begin
@parallel diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz);
update_halo!(T2);
end
T, T2 = T2, T;
end
time_s = toc()
# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(Data.Number); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = time_s/(nt-10); # Execution time per iteration [s]
T_eff = A_eff/t_it; # Effective memory throughput [GB/s]
if (me==0) println("time_s=$time_s T_eff=$T_eff"); end
finalize_global_grid();
end
diffusion3D()Hiding communication with @parallel_indices
The same with a kernel written with @parallel_indices: diffusion3D_multigpucpu_hidecomm_parindices_novis.jl.
Code
julia
const USE_GPU = true
using ImplicitGlobalGrid
using ParallelStencil
@static if USE_GPU
@init_parallel_stencil(CUDA, Float64, 3);
else
@init_parallel_stencil(Threads, Float64, 3);
end
@parallel_indices (ix,iy,iz) function diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz)
if (1<ix<size(T2,1) && 1<iy<size(T2,2) && 1<iz<size(T2,3))
T2[ix,iy,iz] = T[ix,iy,iz] + dt*(Ci[ix,iy,iz]*(
- ((-lam*(T[ix+1,iy,iz] - T[ix,iy,iz])*_dx) - (-lam*(T[ix,iy,iz] - T[ix-1,iy,iz])*_dx))*_dx
- ((-lam*(T[ix,iy+1,iz] - T[ix,iy,iz])*_dy) - (-lam*(T[ix,iy,iz] - T[ix,iy-1,iz])*_dy))*_dy
- ((-lam*(T[ix,iy,iz+1] - T[ix,iy,iz])*_dz) - (-lam*(T[ix,iy,iz] - T[ix,iy,iz-1])*_dz))*_dz)
);
end
return
end
function diffusion3D()
# Physics
lam = 1.0; # Thermal conductivity
c0 = 2.0; # Heat capacity
lx, ly, lz = 1.0, 1.0, 1.0; # Length of computational domain in dimension x, y and z
# Numerics
nx, ny, nz = 512, 512, 512; # Number of gridpoints in dimensions x, y and z
nt = 100; # Number of time steps
me, dims = init_global_grid(nx, ny, nz);
dx = lx/(nx_g()-1); # Space step in x-dimension
dy = ly/(ny_g()-1); # Space step in y-dimension
dz = lz/(nz_g()-1); # Space step in z-dimension
_dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz;
# Array initializations
T = @zeros(nx, ny, nz);
T2 = @zeros(nx, ny, nz);
Ci = @zeros(nx, ny, nz);
# Initial conditions
Ci .= 1/c0; # 1/Heat capacity
T .= 1.7;
T2 .= T; # Assign also T2 to get correct boundary conditions.
# Time loop
dt = min(dx^2,dy^2,dz^2)/lam/maximum(Ci)/8.1; # Time step for 3D Heat diffusion
for it = 1:nt
if (it == 11) tic(); end # Start measuring time.
@hide_communication (16, 2, 2) begin
@parallel diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz);
update_halo!(T2);
end
T, T2 = T2, T;
end
time_s = toc()
# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(Data.Number); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = time_s/(nt-10); # Execution time per iteration [s]
T_eff = A_eff/t_it; # Effective memory throughput [GB/s]
if (me==0) println("time_s=$time_s T_eff=$T_eff"); end
finalize_global_grid();
end
diffusion3D()Hiding communication with boundary conditions
Boundary conditions are set inside the @hide_communication block: diffusion3Dbc_multigpucpu_hidecomm_novis.jl.
Code
julia
const USE_GPU = true
using ImplicitGlobalGrid
import MPI
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
@init_parallel_stencil(CUDA, Float64, 3);
else
@init_parallel_stencil(Threads, Float64, 3);
end
@parallel function diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz)
@inn(T2) = @inn(T) + dt*(lam*@inn(Ci)*(@d2_xi(T)*_dx^2 + @d2_yi(T)*_dy^2 + @d2_zi(T)*_dz^2));
return
end
@parallel function diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz)
@inn(T2) = @inn(T) + dt*(lam*@inn(Ci)*(@d2_xi(T)*_dx^2 + @d2_yi(T)*_dy^2 + @d2_zi(T)*_dz^2));
return
end
function diffusion3D()
# Physics
lam = 1.0; # Thermal conductivity
c0 = 2.0; # Heat capacity
lx, ly, lz = 1.0, 1.0, 1.0; # Length of computational domain in dimension x, y and z
# Numerics
nx, ny, nz = 512, 512, 512; # Number of gridpoints in dimensions x, y and z
nt = 100; # Number of time steps
me, dims = init_global_grid(nx, ny, nz);
dx = lx/(nx_g()-1); # Space step in x-dimension
dy = ly/(ny_g()-1); # Space step in y-dimension
dz = lz/(nz_g()-1); # Space step in z-dimension
_dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz;
# Array initializations
T = @zeros(nx, ny, nz);
T2 = @zeros(nx, ny, nz);
Ci = @zeros(nx, ny, nz);
# Initial conditions
Ci .= 1/c0; # 1/Heat capacity
T .= 1.7;
T2 .= T; # Assign also T2 to get correct boundary conditions.
# Time loop
dt = min(dx^2,dy^2,dz^2)/lam/maximum(Ci)/8.1; # Time step for 3D Heat diffusion
for it = 1:nt
if (it == 11) tic(); end # Start measuring time.
@hide_communication (16, 2, 2) begin
@parallel diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz);
T2[ 1,:,:] .= 1.0;
T2[end,:,:] .= 1.0;
update_halo!(T2);
end
T, T2 = T2, T;
end
time_s = toc()
# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(Data.Number); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = time_s/(nt-10); # Execution time per iteration [s]
T_eff = A_eff/t_it; # Effective memory throughput [GB/s]
if (me==0) println("time_s=$time_s T_eff=$T_eff"); end
finalize_global_grid();
end
diffusion3D()Hiding communication with visualization
The communication is hidden behind computation, and the results are gathered on one process (gather! of ImplicitGlobalGrid.jl) to create the movie shown in the 50-lines multi-xPU example with Plots.jl: diffusion3D_multigpucpu_hidecomm.jl.
Code
julia
const USE_GPU = true
using ImplicitGlobalGrid, Plots
import MPI
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
@init_parallel_stencil(CUDA, Float64, 3);
else
@init_parallel_stencil(Threads, Float64, 3);
end
@parallel function diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz)
@inn(T2) = @inn(T) + dt*(lam*@inn(Ci)*(@d2_xi(T)*_dx^2 + @d2_yi(T)*_dy^2 + @d2_zi(T)*_dz^2));
return
end
function diffusion3D()
# Physics
lam = 1.0; # Thermal conductivity
cp_min = 1.0; # Minimal heat capacity
lx, ly, lz = 10.0, 10.0, 10.0; # Length of computational domain in dimension x, y and z
# Numerics
nx, ny, nz = 256, 256, 256; # Number of gridpoints in dimensions x, y and z
nt = 100000; # Number of time steps
me, dims = init_global_grid(nx, ny, nz);
dx = lx/(nx_g()-1); # Space step in dimension x
dy = ly/(ny_g()-1); # ... in dimension y
dz = lz/(nz_g()-1); # ... in dimension z
_dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz;
# Array initializations
T = @zeros(nx, ny, nz);
T2 = @zeros(nx, ny, nz);
Ci = @zeros(nx, ny, nz);
# Initial conditions (heat capacity and temperature with two Gaussian anomalies each)
Ci .= 1.0./( cp_min .+ Data.Array([5*exp(-((x_g(ix,dx,Ci)-lx/1.5))^2-((y_g(iy,dy,Ci)-ly/2))^2-((z_g(iz,dz,Ci)-lz/1.5))^2) +
5*exp(-((x_g(ix,dx,Ci)-lx/3.0))^2-((y_g(iy,dy,Ci)-ly/2))^2-((z_g(iz,dz,Ci)-lz/1.5))^2) for ix=1:size(T,1), iy=1:size(T,2), iz=1:size(T,3)]) )
T .= Data.Array([100*exp(-((x_g(ix,dx,T)-lx/2)/2)^2-((y_g(iy,dy,T)-ly/2)/2)^2-((z_g(iz,dz,T)-lz/3.0)/2)^2) +
50*exp(-((x_g(ix,dx,T)-lx/2)/2)^2-((y_g(iy,dy,T)-ly/2)/2)^2-((z_g(iz,dz,T)-lz/1.5)/2)^2) for ix=1:size(T,1), iy=1:size(T,2), iz=1:size(T,3)])
T2 .= T; # Assign also T2 to get correct boundary conditions.
# Preparation of visualisation
gr()
ENV["GKSwstype"]="nul"
anim = Animation();
nx_v = (nx-2)*dims[1];
ny_v = (ny-2)*dims[2];
nz_v = (nz-2)*dims[3];
T_v = zeros(nx_v, ny_v, nz_v);
T_nohalo = zeros(nx-2, ny-2, nz-2);
# Time loop
dt = min(dx^2,dy^2,dz^2)*cp_min/lam/8.1; # Time step for the 3D Heat diffusion
for it = 1:nt
if (it == 11) tic(); end # Start measuring time.
if mod(it, 1000) == 1 # Visualize only every 1000th time step
T_nohalo .= Array(T[2:end-1,2:end-1,2:end-1]); # Copy data to CPU removing the halo.
gather!(T_nohalo, T_v) # Gather data on process 0 (could be interpolated/sampled first)
if (me==0) heatmap(transpose(T_v[:,ny_v÷2,:]), aspect_ratio=1); frame(anim); end # Visualize it on process 0.
end
@hide_communication (16, 2, 2) begin
@parallel diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz);
update_halo!(T2);
end
T, T2 = T2, T;
end
time_s = toc()
# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(Data.Number); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = time_s/(nt-10); # Execution time per iteration [s]
T_eff = A_eff/t_it; # Effective memory throughput [GB/s]
if (me==0) println("time_s=$time_s T_eff=$T_eff"); end
# Postprocessing
if (me==0) gif(anim, "diffusion3D.gif", fps = 15) end # Create a gif movie on process 0.
if (me==0) mp4(anim, "diffusion3D.mp4", fps = 15) end # Create a mp4 movie on process 0.
finalize_global_grid(); # Finalize the implicit global grid
end
diffusion3D()Number type as parameter
@init_parallel_stencil is called without numbertype, and the kernels and functions are written generically for any number type (PSNumber): parametrized/diffusion3D_multigpucpu_novis.jl.
Code
julia
const USE_GPU = true
using ImplicitGlobalGrid
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
@init_parallel_stencil(package=CUDA, ndims=3);
else
@init_parallel_stencil(package=Threads, ndims=3);
end
@parallel function diffusion3D_step!(T2::Data.Array{_T}, T::Data.Array{_T}, Ci::Data.Array{_T}, lam::_T, dt::_T, _dx::_T, _dy::_T, _dz::_T) where _T <: PSNumber
@inn(T2) = @inn(T) + dt*(lam*@inn(Ci)*(@d2_xi(T)*_dx^2 + @d2_yi(T)*_dy^2 + @d2_zi(T)*_dz^2));
return
end
function diffusion3D(lam::_T, c0::_T, lx::_T, ly::_T, lz::_T) where _T <: PSNumber
# Numerics
nx, ny, nz = 512, 512, 512; # Number of gridpoints in dimensions x, y and z
nt = 100; # Number of time steps
me, dims = init_global_grid(nx, ny, nz);
dx = lx/(nx_g()-1); # Space step in x-dimension
dy = ly/(ny_g()-1); # Space step in y-dimension
dz = lz/(nz_g()-1); # Space step in z-dimension
_dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz;
# Array initializations
T = @zeros(_T, nx, ny, nz);
T2 = @zeros(_T, nx, ny, nz);
Ci = @zeros(_T, nx, ny, nz);
# Initial conditions
Ci .= 1/c0; # 1/Heat capacity
T .= 1.7;
T2 .= T; # Assign also T2 to get correct boundary conditions.
# Time loop
dt = min(dx^2,dy^2,dz^2)/lam/maximum(Ci)/8.1; # Time step for 3D Heat diffusion
lam, dt, _dx, _dy, _dz = _T.((lam, dt, _dx, _dy, _dz)) # Convert all scalars to the required numbertype.
for it = 1:nt
if (it == 11) tic(); end # Start measuring time.
@parallel diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz);
update_halo!(T2);
T, T2 = T2, T;
end
time_s = toc()
# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(_T); # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = time_s/(nt-10); # Execution time per iteration [s]
T_eff = A_eff/t_it; # Effective memory throughput [GB/s]
if (me==0) println("time_s=$time_s T_eff=$T_eff"); end
finalize_global_grid();
end
################################################################################
# CALL OF THE FUNCTION
# Physics
lam = 1.0; # Thermal conductivity
c0 = 2.0; # Heat capacity
lx, ly, lz = 1.0, 1.0, 1.0; # Length of computational domain in dimension x, y and z
lam, c0, lx, ly, lz = Float64.((lam, c0, lx, ly, lz)) # Convert all scalars to the desired numbertype.
diffusion3D(lam, c0, lx, ly, lz)