How to use data as initial conditions while running on GPU

[iuryt]

Hi,

I couldn't find an example from the docs that uses a regular array, converts to CudaArray and them use it as an initial condition for a GPU run. Do you guys have any example of that that could share? Maybe would be good to consider putting an example like this on the docs.

[glwagner]

Hmm, I believe set!(field, array) should work no matter the architecture of field and the type of array. You can test that with the two-dimensional turbulence example:

https://github.com/CliMA/Oceananigans.jl/blob/main/examples/two_dimensional_turbulence.jl

Raise an issue if that doesn't work!

[iuryt]

It is not working for my example, but I haven't posted because I need a simpler example to reproduce the error. I will use this example and see if it works.

[iuryt]

It worked for the two dimensional turbulence example. So I was checking and comparing the arrays from the example and the one I am trying to use and I found out why I couldn't directly assign the initial conditions while using GPU, but I am still having an error while running the simulation.

For my example, I am reading a netcdf file and it was returning an array with a different format (Noticed it was Array{Float64, .. from the example):

using NCDatasets
ds = Dataset("data/input_coarsened.nc")

ds["b"][:]
> 10×50×48 Array{Union{Missing, Float32}, 3}:

So I simply did:

bᵢ = Float64.(ds["b"][:])

set!(model; b=bᵢ)

and it worked. But, when I try to run the simulation it returns the following error (doesn't happen for CPU):

┌ Info: Initializing simulation...
└ @ Oceananigans.Simulations /home/isimoesdesousa/.julia/packages/Oceananigans/DpyWt/src/Simulations/run.jl:167
┌ Info: i: 0000, t: 0 seconds, Δt: 1.100 minutes, umax = (0.0e+00, 0.0e+00, 0.0e+00) ms⁻¹, wall time: 0 seconds
└ @ Main In[1]:89
┌ Info:     ... simulation initialization complete (38.447 seconds)
└ @ Oceananigans.Simulations /home/isimoesdesousa/.julia/packages/Oceananigans/DpyWt/src/Simulations/run.jl:202
┌ Info: Executing initial time step...
└ @ Oceananigans.Simulations /home/isimoesdesousa/.julia/packages/Oceananigans/DpyWt/src/Simulations/run.jl:112
InvalidIRError: compiling kernel gpu_calculate_Gu!(Cassette.Context{nametype(CUDACtx), Nothing, Nothing, KernelAbstractions.var"##PassType#257", Nothing, Cassette.DisableHooks}, typeof(Oceananigans.Models.NonhydrostaticModels.gpu_calculate_Gu!), KernelAbstractions.CompilerMetadata{KernelAbstractions.NDIteration.StaticSize{(10, 50, 48)}, KernelAbstractions.NDIteration.DynamicCheck, Nothing, Nothing, KernelAbstractions.NDIteration.NDRange{3, KernelAbstractions.NDIteration.StaticSize{(1, 4, 48)}, KernelAbstractions.NDIteration.StaticSize{(16, 16, 1)}, Nothing, Nothing}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, RectilinearGrid{Float64, Periodic, Bounded, Bounded, Float64, Float64, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, Nothing}, UpwindBiasedFifthOrder, FPlane{Float64}, Nothing, Tuple{ScalarDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Oceananigans.TurbulenceClosures.ThreeDimensionalFormulation, typeof(νz), NamedTuple{(:b,), Tuple{typeof(νz)}}}, HorizontalScalarDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, typeof(νh), NamedTuple{(:b,), Tuple{typeof(νh)}}}}, Buoyancy{BuoyancyTracer, Oceananigans.Grids.ZDirection}, NamedTuple{(:velocities, :tracers), Tuple{NamedTuple{(:u, :v, :w), Tuple{Oceananigans.Fields.ZeroField{Int64, 3}, Oceananigans.Fields.ZeroField{Int64, 3}, Oceananigans.Fields.ZeroField{Int64, 3}}}, NamedTuple{(:b,), Tuple{Oceananigans.Fields.ZeroField{Int64, 3}}}}}, NamedTuple{(:u, :v, :w), Tuple{OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}}}, NamedTuple{(:b,), Tuple{OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}}}, Tuple{Nothing, Nothing}, NamedTuple{(:u, :v, :w, :b), NTuple{4, typeof(Oceananigans.Forcings.zeroforcing)}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, NamedTuple{(:time, :iteration, :stage), Tuple{Float64, Int64, Int64}}) resulted in invalid LLVM IR
Reason: unsupported dynamic function invocation (call to overdub)
Stacktrace:
  [1] νz(::Float64, ::Float64, ::Float64, ::Float64)
    @ ./In[1]:35
  [2] overdub
    @ ./In[1]:35
  [3] νᶜᶜᶜ(::Int64, ::Int64, ::Int64, ::RectilinearGrid{Float64, Periodic, Bounded, Bounded, Float64, Float64, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, Nothing}, ::NamedTuple{(:time, :iteration, :stage), Tuple{Float64, Int64, Int64}}, ::Tuple{Center, Center, Center}, ::typeof(νz))
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/closure_kernel_operators.jl:107
  [4] overdub
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/closure_kernel_operators.jl:107
  [5] νᶜᶜᶜ(::Int64, ::Int64, ::Int64, ::RectilinearGrid{Float64, Periodic, Bounded, Bounded, Float64, Float64, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, Nothing}, ::ScalarDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Oceananigans.TurbulenceClosures.ThreeDimensionalFormulation, typeof(νz), NamedTuple{(:b,), Tuple{typeof(νz)}}}, ::Nothing, ::NamedTuple{(:time, :iteration, :stage), Tuple{Float64, Int64, Int64}})
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl:63
  [6] overdub
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl:63
  [7] ν_σᶜᶜᶜ(::Int64, ::Int64, ::Int64, ::RectilinearGrid{Float64, Periodic, Bounded, Bounded, Float64, Float64, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, Nothing}, ::ScalarDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Oceananigans.TurbulenceClosures.ThreeDimensionalFormulation, typeof(νz), NamedTuple{(:b,), Tuple{typeof(νz)}}}, ::Nothing, ::NamedTuple{(:time, :iteration, :stage), Tuple{Float64, Int64, Int64}}, ::typeof(Oceananigans.TurbulenceClosures.Σ₁₁), ::OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, ::OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, ::OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}})
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/closure_kernel_operators.jl:54
  [8] overdub
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/closure_kernel_operators.jl:54
  [9] viscous_flux_ux(::Int64, ::Int64, ::Int64, ::RectilinearGrid{Float64, Periodic, Bounded, Bounded, Float64, Float64, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, Nothing}, ::ScalarDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Oceananigans.TurbulenceClosures.ThreeDimensionalFormulation, typeof(νz), NamedTuple{(:b,), Tuple{typeof(νz)}}}, ::Nothing, ::NamedTuple{(:u, :v, :w), Tuple{OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}}}, ::NamedTuple{(:b,), Tuple{OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}}}, ::NamedTuple{(:time, :iteration, :stage), Tuple{Float64, Int64, Int64}}, ::Buoyancy{BuoyancyTracer, Oceananigans.Grids.ZDirection})
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl:89
 [10] overdub
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl:89
 [11] viscous_flux_ux(::Int64, ::Int64, ::Int64, ::RectilinearGrid{Float64, Periodic, Bounded, Bounded, Float64, Float64, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, Nothing}, ::Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, ::ScalarDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Oceananigans.TurbulenceClosures.ThreeDimensionalFormulation, typeof(νz), NamedTuple{(:b,), Tuple{typeof(νz)}}}, ::Nothing, ::NamedTuple{(:u, :v, :w), Tuple{OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}}}, ::NamedTuple{(:b,), Tuple{OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}}}, ::NamedTuple{(:time, :iteration, :stage), Tuple{Float64, Int64, Int64}}, ::Buoyancy{BuoyancyTracer, Oceananigans.Grids.ZDirection})
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/implicit_explicit_time_discretization.jl:46
 [12] overdub
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/implicit_explicit_time_discretization.jl:46
 [13] _viscous_flux_ux(::Int64, ::Int64, ::Int64, ::RectilinearGrid{Float64, Periodic, Bounded, Bounded, Float64, Float64, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, Nothing}, ::Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, ::ScalarDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Oceananigans.TurbulenceClosures.ThreeDimensionalFormulation, typeof(νz), NamedTuple{(:b,), Tuple{typeof(νz)}}}, ::Nothing, ::NamedTuple{(:u, :v, :w), Tuple{OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}}}, ::NamedTuple{(:b,), Tuple{OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}}}, ::NamedTuple{(:time, :iteration, :stage), Tuple{Float64, Int64, Int64}}, ::Buoyancy{BuoyancyTracer, Oceananigans.Grids.ZDirection})
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/closure_kernel_operators.jl:4
 [14] overdub
    @ ~/.julia/packages/Oceananigans/DpyWt/src/TurbulenceClosures/closure_kernel_operators.jl:4
 [15] Ax_qᶜᶜᶜ(::Int64, ::Int64, ::Int64, ::RectilinearGrid{Float64, Periodic, Bounded, Bounded, Float64, Float64, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, CUDA.CuDeviceVector{Float64, 1}}, Nothing}, ::typeof(Oceananigans.TurbulenceClosures._viscous_flux_ux), ::Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, ::ScalarDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Oceananigans.TurbulenceClosures.ThreeDimensionalFormulation, typeof(νz), NamedTuple{(:b,), Tuple{typeof(νz)}}}, ::Nothing, ::NamedTuple{(:u, :v, :w), Tuple{OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}, OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}}}, ::NamedTuple{(:b,), Tuple{OffsetArrays.OffsetArray{Float64, 3, CUDA.CuDeviceArray{Float64, 3, 1}}}}, ::NamedTuple{(:time, :iteration, :stage), Tuple{Float64, Int64, Int64}}, ::Buoyancy{BuoyancyTracer, Oceananigans.Grids.ZDirection})
    @ ~/.julia/packages/Oceananigans/DpyWt/src/Operators/products_between_fields_and_grid_metrics.jl:12
.
.
.
.

Any idea on what is happening or should I try to share the full code with the NC file?
Anyway I will be exploring this tomorrow to see if I can find anything else.

[glwagner]

It looks like the problem doesn't have to do with the initial condition, but instead with the function νz(::Float64, ::Float64, ::Float64, ::Float64). Does that function reference any object as a global variable? You may be missing const somewhere.

[iuryt]

You are right!

The problem was here, I didn't have const before

using NCDatasets
ds = Dataset("data/input_coarsened.nc")
const Nx, Ny, Nz = size(ds["b"])