KernelFunctionOperation question

[SamDeAbreu]

Hi everyone,

I'm trying to compute a quantity that involves a lot of derivatives of the velocity fields on GPU over a large grid. The quantity is:
Q = 0.5*((u_y - v_x)**2 + (u_z - w_x)**2 + (v_z - w_y)**2 - (u_x**2 + v_y**2 + w_z**2) - (u_y + v_x)**2 - (u_z + w_x)**2 - (v_z+w_y)**2)
evaluated at the surface. I originally computed each derivative separately, used @at (Center, Center, Center) at each one, and then computed Q. This resulted in a segmentation fault, which I believe was a result of running out of memory from storing all of the derivative fields.

I decided then that maybe it is time to learn about kernel operations to avoid storing all the fields. This is how I went about it:

u = model.velocities.u
v = model.velocities.v
w = model.velocities.w
function Q_ccc(i, j, k, grid, u, v, w)
    u_x = ∂xᶜᵃᵃ(i, j, k, grid, u) # F, C, C -> C, C, C
    u_y = ℑxyᶜᶜᵃ(i, j, k, grid, ∂yᵃᶠᵃ, u) # F, C, C -> F, F, C -> C, C, C
    u_z = ℑxzᶜᵃᶜ(i, j, k, grid, ∂zᵃᵃᶠ, u) # F, C, C -> F, C, F -> C, C, C

    v_x = ℑxyᶜᶜᵃ(i, j, k, grid, ∂xᶠᵃᵃ, v) # C, F, C -> F, F, C -> C, C, C
    v_y = ∂yᵃᶜᵃ(i, j, k, grid, v) # C, F, C -> C, C, C
    v_z = ℑyzᵃᶜᶜ(i, j, k, grid, ∂zᵃᵃᶠ, v) # C, F, C -> C, F, F -> C, C, C 

    w_x = ℑxzᶜᵃᶜ(i, j, k, grid, ∂xᶠᵃᵃ, w) # C, C, F -> F, C, F -> C, C, C
    w_y = ℑyzᵃᶜᶜ(i, j, k, grid, ∂yᵃᶠᵃ, w) # C, C, F -> C, F, F -> C, C, C
    w_z = ∂zᵃᵃᶜ(i, j, k, grid, w)
    return 0.5/U_b^2*((u_y - v_x)^2 + (u_z - w_x)^2 + (v_z - w_y)^2 - (u_x^2 + v_y^2 + w_z^2) - (u_y + v_x)^2 - (u_z + w_x)^2 - (v_z + w_y)^2)
end
Q_op = KernelFunctionOperation{Center, Center, Center}(Q_ccc, grid, u, v, w)
Q = Field(Q_op)

Is this correct?

[glwagner]

Without inspecting very closely to find a typo, the methodology is correct. However, it is "more correct" to delay averaging as long as possible. For example, your expression has a pattern:

Q1_ffc = (u_y - v_x)^2 - (u_y + v_x)^2
Q2_fcf = (u_z - w_x)^2 - (u_z + w_x)^2
Q3_cff = (v_z - w_y)^2 - (v_z + w_y)^2
Q4_ccc = (u_x^2 + v_y^2 + w_z^2)

You could compute all of these in separate AbstractOperation, and then finally assemble into Q. You can also use a function, but in that case I would use this pattern:

Q1_ffc(i, j, k, grid, u, v) = (∂yᵃᶠᵃ(i, j, k, grid, u) - ∂xᶠᵃᵃ(i, j, k, grid, v))^2 - (∂yᵃᶠᵃ(i, j, k, grid, u) + ∂xᶠᵃᵃ(i, j, k, grid, v))^2

Q_ccc(i, j, k, grid, u, v, w) = ℑxyᶜᶜᵃ(i, j, k, grid, Q1_ffc, u, v) + # etc

in other words, form the components at the locations they "naturally" occur at and then average at the end.

You also shouldn't get a segmentation fault. Is U_b a constant?

[glwagner]

Does this work?

Q = @at (Center, Center, Center) Field(0.5*(term1 + term3)/U_b^2)

I think you need to apply @at to the operation, eg

Q_op = @at (Center, Center, Center) (term1 + term3) / 2Ub^2
Q = Field(Q_op)

[SamDeAbreu]

That seems to produce consistent seg faults

[glwagner]

Ok new idea. Can you try using a larger halo for the grid? Perhaps halo = (7, 7, 7) to be conservative, but halo = (4, 4, 4) might simply work (the default is halo = (3, 3, 3)). Yes, good to use the version of the code that produces a consistent error 😂

[SamDeAbreu]

That seems to have fixed it? No more seg faults now. Thank you!

[glwagner]

Idea for the future might be to figure out how to check how big of a halo an operation will need. I think this can be done (from a logic perspective) though it may be a little intricate (eg we have to compute the net halo region size using interpolations+differences I think) cc @simone-silvestri