How to use ∂ function in my code or how to compute the mixed second-order partial derivatives or first-order cross partial derivatives of the velocity field **u**?

’‘’
using WaterLily,StaticArrays,GLMakie,Makie,Plots
using LinearAlgebra: eigvals, Hermitian

function donut(L;Re=1e3,mem=Array,U=1)
center,r = SA[3L,3L,0], L/2
bodya = AutoBody() do xyz,t
x,y,z = xyz - center
√sum(abs2,SA[x,y,0])-r
end
cycenter=SA[3L,3L,4L]
bodyb = AutoBody() do xyz,t
x,y,z = xyz - cycenter
a=√sum(abs2,SA[x,y,0])-L
b=√sum(abs2,SA[0,0,z])-r
max(a,b)
end
body=bodya+bodyb
Simulation((8L,6L,8L),(U,0,0),L;ν=U*L/Re,body,mem)
end

function Q(I::CartesianIndex{3}, u)
J = @SMatrix [∂(i, j, I, u) for i in 1:3, j in 1:3]
S = (J + J’) / 2
Ω = (J - J’) / 2
norm_S_sq = sum(S.^2)
norm_Ω_sq = sum(Ω.^2)

0.5 * (norm_Ω_sq - norm_S_sq)                        

end

function ω!(arr, sim)
a = sim.flow.σ
WaterLily.@inside a[I] = Q(I,sim.flow.u)
copyto!(arr, a[inside(a)]) # copy to CPU
end

sim = donut(2^3);
t₀ = sim_time(sim)
duration = 100.0
step = 0.25
viz!(sim;f=ω!,duration,step,framerate=5,algorithm=:mip,colormap=:algae)
‘’‘

The code above is what I want to achieve. I saw the following code in the StaticArrays package.
‘’’
function λ₂(I::CartesianIndex{3},u)
J = @SMatrix [∂(i,j,I,u) for i ∈ 1:3, j ∈ 1:3]
S,Ω = (J+J’)/2,(J-J’)/2
eigvals(Hermitian(S^2+Ω^2))[2]
end
‘’’
So I want to write a function following this example. But when I run my code, an error occurs. The error is as follows.
‘’’
1-element ExceptionStack:
LoadError: TaskFailedException

nested task error: UndefVarError: `∂` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Stacktrace:
 [1] (::var"#Q##2#Q##3"{CartesianIndex{3}, Array{Float32, 4}})(i::Int64, j::Int64)
   @ Main D:\julia\julia_pakage\packages\StaticArrays\DsPgf\src\SMatrix.jl:63
 [2] macro expansion
   @ D:\julia\julia_pakage\packages\StaticArrays\DsPgf\src\SMatrix.jl:64 [inlined]
 [3] Q(I::CartesianIndex{3}, u::Array{Float32, 4})
   @ Main f:\julia_code\3D阶梯圆柱\threed.jl:36
 [4] macro expansion
   @ D:\julia\julia_pakage\packages\WaterLily\O8g3K\src\util.jl:146 [inlined]
 [5] macro expansion
   @ D:\julia\julia_pakage\packages\KernelAbstractions\X5fk1\src\macros.jl:314 [inlined]
 [6] (::var"#cpu_##kern_#286#ω!##5")(__ctx__::KernelAbstractions.CompilerMetadata{KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIteration.NoDynamicCheck, CartesianIndex{3}, CartesianIndices{3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}, KernelAbstractions.NDIteration.NDRange{3, KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIteration.StaticSize{(64, 1, 1)}, CartesianIndices{3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}, Nothing}}, a::Array{Float32, 3}, u::Array{Float32, 4}, I0::CartesianIndex{3})
   @ Main .\none:0
 [7] __thread_run(tid::Int64, len::Int64, rem::Int64, obj::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(64,)}, KernelAbstractions.NDIteration.DynamicSize, var"#cpu_##kern_#286#ω!##5"}, ndrange::Tuple{Int64, Int64, Int64}, iterspace::KernelAbstractions.NDIteration.NDRange{3, KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIteration.StaticSize{(64, 1, 1)}, CartesianIndices{3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}, Nothing}, args::Tuple{Array{Float32, 3}, Array{Float32, 4}, CartesianIndex{3}}, dynamic::KernelAbstractions.NDIteration.NoDynamicCheck)
   @ KernelAbstractions D:\julia\julia_pakage\packages\KernelAbstractions\X5fk1\src\cpu.jl:145
 [8] (::KernelAbstractions.var"#__run##4#__run##5"{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(64,)}, KernelAbstractions.NDIteration.DynamicSize, var"#cpu_##kern_#286#ω!##5"}, Tuple{Int64, Int64, Int64}, KernelAbstractions.NDIteration.NDRange{3, KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIteration.StaticSize{(64, 1, 1)}, CartesianIndices{3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}, Nothing}, Tuple{Array{Float32, 3}, Array{Float32, 4}, CartesianIndex{3}}, KernelAbstractions.NDIteration.NoDynamicCheck, Int64})()
   @ KernelAbstractions D:\julia\julia_pakage\packages\KernelAbstractions\X5fk1\src\cpu.jl:120

…and 19 more exceptions.

Stacktrace:
[1] sync_end(c::Channel{Any})
@ Base .\task.jl:601
[2] macro expansion
@ .\task.jl:634 [inlined]
[3] run(obj::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(64,)}, KernelAbstractions.NDIteration.DynamicSize, var"cpu##kern#286#ω!##5"}, ndrange::Tuple{Int64, Int64, Int64}, iterspace::KernelAbstractions.NDIteration.NDRange{3, KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIteration.StaticSize{(64, 1, 1)}, CartesianIndices{3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}, Nothing}, args::Tuple{Array{Float32, 3}, Array{Float32, 4}, CartesianIndex{3}}, dynamic::KernelAbstractions.NDIteration.NoDynamicCheck, static_threads::Bool)
@ KernelAbstractions D:\julia\julia_pakage\packages\KernelAbstractions\X5fk1\src\cpu.jl:119
[4] (::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(64,)}, KernelAbstractions.NDIteration.DynamicSize, var"#cpu_##kern_#286#ω!##5"})(::Array{Float32, 3}, ::Vararg{Any}; ndrange::Tuple{Int64, Int64, Int64}, workgroupsize::Nothing)
@ KernelAbstractions D:\julia\julia_pakage\packages\KernelAbstractions\X5fk1\src\cpu.jl:46
[5] (::var"##kern#285#ω!##9"{Simulation})(a::Array{Float32, 3}, u::Array{Float32, 4})
@ Main D:\julia\julia_pakage\packages\WaterLily\O8g3K\src\util.jl:149
[6] macro expansion
@ D:\julia\julia_pakage\packages\WaterLily\O8g3K\src\util.jl:151 [inlined]
[7] ω!(arr::Array{Float32, 3}, sim::Simulation)
@ Main f:\julia_code\3D阶梯圆柱\threed.jl:56
[8] viz!(sim::Simulation; f::typeof(ω!), duration::Float64, step::Float64, remeasure::Bool, verbose::Bool, λ::Function, udf::Nothing, udf_kwargs::Nothing, d::Int64, CIs::Nothing, cut::Nothing, tidy_colormap::Bool, body::Bool, body_color::Symbol, body2mesh::Bool, video::Nothing, hidedecorations::Bool, elevation::Float64, azimuth::Float64, framerate::Int64, compression::Int64, theme::Nothing, fig_size::Nothing, fig_pad::Int64, fig::Nothing, ax::Nothing, kwargs::@Kwargs{algorithm::Symbol, colormap::Symbol})
@ WaterLilyMakieExt D:\julia\julia_pakage\packages\WaterLily\O8g3K\ext\WaterLilyMakieExt.jl:184
[9] top-level scope
@ f:\julia_code\3D阶梯圆柱\threed.jl:66
[10] eval(m::Module, e::Any)
@ Core .\boot.jl:489
[11] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
@ Base .\loading.jl:2842
[12] include_string(m::Module, txt::String, fname::String)
@ Base .\loading.jl:2852
[13] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
@ VSCodeServer c:\Users\intel.vscode\extensions\julialang.language-julia-1.158.2\scripts\packages\VSCodeServer\src\eval.jl:271
[14] (::VSCodeServer.var"#repl_runcode_request##6#repl_runcode_request##7"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer c:\Users\intel.vscode\extensions\julialang.language-julia-1.158.2\scripts\packages\VSCodeServer\src\eval.jl:181
[15] withpath(f::VSCodeServer.var"#repl_runcode_request##6#repl_runcode_request##7"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
@ VSCodeServer c:\Users\intel.vscode\extensions\julialang.language-julia-1.158.2\scripts\packages\VSCodeServer\src\repl.jl:278
[16] (::VSCodeServer.var"#repl_runcode_request##4#repl_runcode_request##5"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer c:\Users\intel.vscode\extensions\julialang.language-julia-1.158.2\scripts\packages\VSCodeServer\src\eval.jl:179
[17] hideprompt(f::VSCodeServer.var"#repl_runcode_request##4#repl_runcode_request##5"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
@ VSCodeServer c:\Users\intel.vscode\extensions\julialang.language-julia-1.158.2\scripts\packages\VSCodeServer\src\repl.jl:38
[18] #repl_runcode_request##2
@ c:\Users\intel.vscode\extensions\julialang.language-julia-1.158.2\scripts\packages\VSCodeServer\src\eval.jl:150 [inlined]
[19] with_logstate(f::VSCodeServer.var"#repl_runcode_request##2#repl_runcode_request##3"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, logstate::Base.CoreLogging.LogState)
@ Base.CoreLogging .\logging\logging.jl:540
[20] with_logger
@ .\logging\logging.jl:651 [inlined]
[21] (::VSCodeServer.var"#repl_runcode_request##0#repl_runcode_request##1"{VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer c:\Users\intel.vscode\extensions\julialang.language-julia-1.158.2\scripts\packages\VSCodeServer\src\eval.jl:263
[22] (::VSCodeServer.var"#start_eval_backend##0#start_eval_backend##1")()
@ VSCodeServer c:\Users\intel.vscode\extensions\julialang.language-julia-1.158.2\scripts\packages\VSCodeServer\src\eval.jl:34
in expression starting at f:\julia_code\3D阶梯圆柱\threed.jl:66
‘’’
So I want to konw how can I use ∂ in my code or how can I compute the mixed second-order partial derivatives or first-order cross partial derivatives of the velocity field u, using a central difference scheme.
I am looking forward to your answer. Thank you very much.

Are you doing this as part of a university course or something? The same question was asked yesterday:

1 Like

I use it to visualize three-dimensional vorticity fields, which constitutes an important part of my scientific research.

Where do you expect the \partial function to come from? It isn’t defined anywhere in your code. Usually, these computations are performed with an automatic differentiation library, so unless you have an explicit formula for ∂(i, j, I, u), you may need to pick one of these libraries in the Julia ecosystem.

Also, please try to format your posts using the correct backticks ```, so that the code is displayed properly, see Please read: make it easier to help you.

2 Likes

Thank you for your suggestion. I found ∂(i, j, I, u) in the path D:\julia\julia_pakage\packages\WaterLily\08g3K\src\Metrics.jl, which contains the definition of ∂(i, j, I, u) as follows.

"""
    ∂(i,j,I,u)

Compute ``∂uᵢ/∂xⱼ`` at center of cell `I`. Cross terms are computed
less accurately than inline terms because of the staggered grid.
"""
@fastmath @inline (i,j,I,u) = (i==j ? ∂(i,I,u) :
        @inbounds(u[I+δ(j,I),i]+u[I+δ(j,I)+δ(i,I),i]
                 -u[I-δ(j,I),i]-u[I-δ(j,I)+δ(i,I),i])/4)