Max Köhler (Jun 14 2021 at 08:12): Max Köhler (Jun 14 2021 at 08:12):It's time to get my code type stable and I'm wondering whats going on in this snippet:
if !(convex)
convex = false
relaxed = true
ξ = volumefraction(state, F, F⁺,F¯)
if (0.0 < ξ < 1.0) && (F[1] >= F⁺[1] || F[1] <= F¯[1])
relaxed = false
else
relaxed = true
end
this part gives me the following code_warntype
9 ┄─ %75 = !convex::Bool
└─── goto #39 if not %75
10 ─ (convex = false)
│ (relaxed = true)
│ %79 = ConvexDamage.volumefraction(state, F, F⁺, F¯)::Float64
│ Core.setfield!(ξ@_15, :contents, %79)
│ %81 = Core.isdefined(ξ@_15, :contents)::Bool
└─── goto #12 if not %81
11 ─ goto #13
12 ─ Core.NewvarNode(:(ξ@_28))
└─── ξ@_28
13 ┄ %86 = Core.getfield(ξ@_15, :contents)::Any
│ %87 = (0.0 < %86)::Any
Now I'm wondering, why Core.getfield(ξ@_15, :contents) is of type ::Any if ξ is correctly interferred to a Float64 before?
Kwaku Oskin (Jun 14 2021 at 08:16):Is it a part of a function? Or all of it is defined in global space?
Kwaku Oskin (Jun 14 2021 at 08:18):Maybe somewhere else ξ is changed to a different type.
Sebastian Pfitzner (Jun 14 2021 at 08:19):Also, what's that Core.isdefined doing in there?
Max Köhler (Jun 14 2021 at 08:19):Everything is inside a function
function constitutive_driver(F, material::RelaxedDamage, state::RelaxedDamageState)
state.relaxed ? convex = false : convex = true
W_min = state.W_min
W = state.W
F⁺, F¯ = state.damage⁺.Fₖ, state.damage¯.Fₖ
if convex
W_min, F⁺, F¯ = convexify(F, material, state)
ξ = (F[1] - F¯[1]) / (F⁺[1] - F¯[1])
d = (F⁺[1] - F¯[1]) / F[1]
W = W_energy(F, material, state.damage)
convex = W < W_min || isapprox(W_min,W)
end
if !(convex)
convex = false
relaxed = true
ξ = volumefraction(state, F, F⁺,F¯)
if (0.0 < ξ < 1.0) && (F[1] >= F⁺[1] || F[1] <= F¯[1])
relaxed = false
else
relaxed = true
end
The first ξ assignment is interfered to a Float64 as well
Max Köhler (Jun 14 2021 at 08:20):I think Core.isdefined is related to the upper assignment in the sense of checking if it's already assigned, but I'm not sure.
Max Köhler (Jun 14 2021 at 08:21):this is the code_warntype of the upper part related to xi
│ %47 = Base.getindex(F, 1)::Float64
│ %48 = Base.getindex(F¯, 1)::Float64
│ %49 = (%47 - %48)::Float64
│ %50 = Base.getindex(F⁺, 1)::Float64
│ %51 = Base.getindex(F¯, 1)::Float64
│ %52 = (%50 - %51)::Float64
│ %53 = (%49 / %52)::Float64
│ Core.setfield!(ξ@_15, :contents, %53)
Ah, I see, so xi is getting Core.Boxed
Sebastian Pfitzner (Jun 14 2021 at 08:23):You could try a manual local xi::Float64 (or whatever type it should have) before the first condition
Max Köhler (Jun 14 2021 at 08:23):What is Core.Box ? And more importantly how can I get rid of it :D
Max Köhler (Jun 14 2021 at 08:23):ah ok, let me try
Sebastian Pfitzner (Jun 14 2021 at 08:24):And it might also help to rewrite that to an if else block; I'm not sure if Julia can optimize a if cond/if !cond block the same way
Max Köhler (Jun 14 2021 at 08:33):Did you mean something like this ?
if !(convex)
convex = false
relaxed = true
local ξ::Float64 = volumefraction(state, F, F⁺,F¯)
if (0.0 < ξ < 1.0) && (F[1] >= F⁺[1] || F[1] <= F¯[1])
relaxed = false
else
relaxed = true
end
This is still Core.Box according to @code_warntype
Max Köhler (Jun 14 2021 at 08:47):Outcommenting the first assignment should also fix it, shouldn't it? Unfortunately, it doesn't
function constitutive_driver(F, material::RelaxedDamage, state::RelaxedDamageState)
state.relaxed ? convex = false : convex = true
W_min = state.W_min
W = state.W
F⁺, F¯ = state.damage⁺.Fₖ, state.damage¯.Fₖ
if convex
W_min, F⁺, F¯ = convexify(F, material, state)
#ξ = (F[1] - F¯[1]) / (F⁺[1] - F¯[1])
#d = (F⁺[1] - F¯[1]) / F[1]
W = W_energy(F, material, state.damage)
convex = W < W_min || isapprox(W_min,W)
end
if !(convex)
convex = false
relaxed = true
ξ = volumefraction(state, F, F⁺,F¯)
if (0.0 < ξ < 1.0) && (F[1] >= F⁺[1] || F[1] <= F¯[1])
relaxed = false
else
relaxed = true
end
no, I meant something like
function constitutive_driver(F, material::RelaxedDamage, state::RelaxedDamageState)
state.relaxed ? convex = false : convex = true
W_min = state.W_min
W = state.W
F⁺, F¯ = state.damage⁺.Fₖ, state.damage¯.Fₖ
local ξ::Float64
if convex
W_min, F⁺, F¯ = convexify(F, material, state)
ξ = (F[1] - F¯[1]) / (F⁺[1] - F¯[1])
d = (F⁺[1] - F¯[1]) / F[1]
W = W_energy(F, material, state.damage)
convex = W < W_min || isapprox(W_min,W)
end
if !(convex)
convex = false
relaxed = true
ξ = volumefraction(state, F, F⁺,F¯)
if (0.0 < ξ < 1.0) && (F[1] >= F⁺[1] || F[1] <= F¯[1])
relaxed = false
else
relaxed = true
end
or just ξ = 0.0, I guess
Sebastian Pfitzner (Jun 14 2021 at 10:17):also, disregard that comment about using an if/else block, I didn't read your code correctly the first time
Max Köhler (Jun 14 2021 at 10:39):Doing so:
function constitutive_driver(F, material::RelaxedDamage, state::RelaxedDamageState)
state.relaxed ? convex = false : convex = true
W_min = state.W_min
W = state.W
F⁺, F¯ = state.damage⁺.Fₖ, state.damage¯.Fₖ
local ξ::Float64
local d::Float64
if convex
W_min, F⁺, F¯ = convexify(F, material, state)
ξ = (F[1] - F¯[1]) / (F⁺[1] - F¯[1])
d = (F⁺[1] - F¯[1]) / F[1]
W = W_energy(F, material, state.damage)
convex = W < W_min || isapprox(W_min,W)
end
if !(convex)
convex = false
relaxed = true
ξ = volumefraction(state, F, F⁺,F¯)
if (0.0 < ξ < 1.0) && (F[1] >= F⁺[1] || F[1] <= F¯[1])
relaxed = false
else
relaxed = true
end
fixes pretty much all propagated type instabilities, however the first xi and d statement are now a Core.Box (see %54 and %63)
5 ── %38 = ConvexDamage.convexify(F, material, state)::Tuple{Float64, Tensor{2, 1, Float64, 1}, Tensor{2, 1, Float64, 1}}
│ %39 = Base.indexed_iterate(%38, 1)::Core.PartialStruct(Tuple{Float64, Int64}, Any[Float64, Core.Const(2)])
│ (W_min = Core.getfield(%39, 1))
│ (@_9 = Core.getfield(%39, 2))
│ %42 = Base.indexed_iterate(%38, 2, @_9::Core.Const(2))::Core.PartialStruct(Tuple{Tensor{2, 1, Float64, 1}, Int64}, Any[Tensor{2, 1, Float64, 1}, Core.Const(3)])
│ (F⁺ = Core.getfield(%42, 1))
│ (@_9 = Core.getfield(%42, 2))
│ %45 = Base.indexed_iterate(%38, 3, @_9::Core.Const(3))::Core.PartialStruct(Tuple{Tensor{2, 1, Float64, 1}, Int64}, Any[Tensor{2, 1, Float64, 1}, Core.Const(4)])
│ (F¯ = Core.getfield(%45, 1))
│ %47 = Base.getindex(F, 1)::Float64
│ %48 = Base.getindex(F¯, 1)::Float64
│ %49 = (%47 - %48)::Float64
│ %50 = Base.getindex(F⁺, 1)::Float64
│ %51 = Base.getindex(F¯, 1)::Float64
│ %52 = (%50 - %51)::Float64
│ %53 = (%49 / %52)::Float64
│ %54 = ξ@_11::Core.Box
│ %55 = Base.convert(ConvexDamage.Float64, %53)::Float64
│ %56 = Core.typeassert(%55, ConvexDamage.Float64)::Float64
│ Core.setfield!(%54, :contents, %56)
│ %58 = Base.getindex(F⁺, 1)::Float64
│ %59 = Base.getindex(F¯, 1)::Float64
│ %60 = (%58 - %59)::Float64
│ %61 = Base.getindex(F, 1)::Float64
│ %62 = (%60 / %61)::Float64
│ %63 = d@_10::Core.Box
which corresponds to these lines
W_min, F⁺, F¯ = convexify(F, material, state)
ξ = (F[1] - F¯[1]) / (F⁺[1] - F¯[1])
d = (F⁺[1] - F¯[1]) / F[1]
Sorry for asking now several times, but I don't really understand what's going on. I would understand the compiler complains if the operation before would yield Any or something else, but code_warntype explicitely says it's Float64, so I dont' really get it
Sebastian Pfitzner (Jun 14 2021 at 10:58):Did you try initializing xi and d to e.g. 0.0? Also, maybe try a type annotation for the re-assignments
Max Köhler (Jun 14 2021 at 11:02):If I try the initialization with 0.0 i get already at initialization Core.Box
│ %24 = Base.getproperty(state, :relaxed)::Bool
└─── goto #3 if not %24
2 ── (convex = false)
└─── goto #4
3 ── (convex = true)
4 ┄─ (W_min = Base.getproperty(state, :W_min))
│ (W = Base.getproperty(state, :W))
│ %31 = Base.getproperty(state, :damage⁺)::ConvexDamage.Damage{1, Float64, 1}
│ %32 = Base.getproperty(%31, :Fₖ)::Tensor{2, 1, Float64, 1}
│ %33 = Base.getproperty(state, :damage¯)::ConvexDamage.Damage{1, Float64, 1}
│ %34 = Base.getproperty(%33, :Fₖ)::Tensor{2, 1, Float64, 1}
│ (F⁺ = %32)
│ (F¯ = %34)
│ %37 = ξ@_15::Core.Box
│ %38 = Base.convert(ConvexDamage.Float64, 0.0)::Core.Const(0.0)
│ %39 = Core.typeassert(%38, ConvexDamage.Float64)::Core.Const(0.0)
│ Core.setfield!(%37, :contents, %39)
│ %41 = d@_14::Core.Box
corresponds to the code part above the first if statement
Kwaku Oskin (Jun 14 2021 at 11:02):Just out of curiousity, if you use ξ2 instead of ξ (and comment out everything after conditions), will it still be a Box?
I mean, something like
function constitutive_driver(F, material::RelaxedDamage, state::RelaxedDamageState)
state.relaxed ? convex = false : convex = true
W_min = state.W_min
W = state.W
F⁺, F¯ = state.damage⁺.Fₖ, state.damage¯.Fₖ
local ξ::Float64
local d::Float64
if convex
W_min, F⁺, F¯ = convexify(F, material, state)
ξ = (F[1] - F¯[1]) / (F⁺[1] - F¯[1])
d = (F⁺[1] - F¯[1]) / F[1]
W = W_energy(F, material, state.damage)
convex = W < W_min || isapprox(W_min,W)
end
if !(convex)
convex = false
relaxed = true
ξ2 = volumefraction(state, F, F⁺,F¯)
if (0.0 < ξ2 < 1.0) && (F[1] >= F⁺[1] || F[1] <= F¯[1])
relaxed = false
else
relaxed = true
end
end
end
Note, that function ends here (there is no continuation, which can introduce additional difficulties for compiler).
Max Köhler (Jun 14 2021 at 11:04):In general this is the full function, if it matters:
function constitutive_driver(F, material::RelaxedDamage, state::RelaxedDamageState)
state.relaxed ? convex = false : convex = true
W_min = state.W_min
W = state.W
F⁺, F¯ = state.damage⁺.Fₖ, state.damage¯.Fₖ
ξ = 0.0
d = 0.0
if convex
W_min, F⁺, F¯ = convexify(F, material, state)
ξ::Float64 = (F[1] - F¯[1]) / (F⁺[1] - F¯[1])
d::Float64 = (F⁺[1] - F¯[1]) / F[1]
W = W_energy(F, material, state.damage)
convex = W < W_min || isapprox(W_min,W)
end
if !(convex)
relaxed = true
ξ = volumefraction(state, F, F⁺,F¯)
if (0.0 < ξ < 1.0) && (F[1] >= F⁺[1] || F[1] <= F¯[1])
relaxed = false
else
relaxed = true
end
d = ξ != 1 ? (F⁺[1] - F¯[1]) / F[1] : 0.0
#P, 𝔸, ψ⁺, ψ¯ = homogenized_derivatives(ξ, d, F, material, state) FIXME
𝔸, P = Tensors.hessian(
y -> W_homogenized(ξ, d, y, material, state),
F,
:all)
ψ⁺ = Ψ(tdot(F⁺), material.base_material)
ψ¯ = Ψ(tdot(F¯), material.base_material)
σ = symmetric(inv(det(F)) * P ⋅ F')
W = W_homogenized(ξ, d, F, material, state)
unrelaxed = Damage(ψₖ = state.damage.ψₖ, βₖ = state.damage.βₖ, Fₖ=F)
phase⁺ = Damage(ψₖ = ψ⁺, βₖ = max(state.damage⁺.βₖ, ψ⁺), Fₖ=F⁺)
phase¯ = Damage(ψₖ = ψ¯, βₖ = max(state.damage¯.βₖ, ψ¯), Fₖ=F¯)
else
#𝔸, P, ψ = strain_energy_derivatives(F, material, state)
𝔸, P = Tensors.hessian(y -> W_energy(y, material, state.damage), F, :all)
ψ = Ψ(tdot(F), material.base_material)
relaxed = false
σ = symmetric(inv(det(F)) * state.P ⋅ F')
ξ = state.ξₖ
d = state.dₖ
unrelaxed = Damage(ψₖ = ψ, βₖ = max(state.damage.βₖ, ψ), Fₖ=F)
phase⁺ = Damage(ψₖ = ψ, βₖ = max(state.damage.βₖ, ψ))
phase¯ = Damage(ψₖ = ψ, βₖ = max(state.damage.βₖ, ψ))
end
newstate = RelaxedDamageState(unrelaxed,phase⁺,phase¯,P,σ, ξ,d, W_min, W, convex, relaxed)
return P, 𝔸, newstate
end
I'll try it with a second xi variable
Sebastian Pfitzner (Jun 14 2021 at 11:06):I feel like you'd need to annotate all xi assignments if you want to get rid of the box
Sebastian Pfitzner (Jun 14 2021 at 11:06):but my intuition for that isn't great either apart from closures :stuck_out_tongue:
Kwaku Oskin (Jun 14 2021 at 11:08):Oh, there is another ξ in ξ = state.ξₖ.
Kwaku Oskin (Jun 14 2021 at 11:08):Probably this is where you get type instability.
Max Köhler (Jun 14 2021 at 11:09):If I annotate multiple ones I get multiple type annotation error, am I doing it wrong by just writing ξ::Float64 = ... in front of each assignment?
Sebastian Pfitzner (Jun 14 2021 at 11:09):ah, yeah, you'd want to annotate the rhs, not the lhs
Max Köhler (Jun 14 2021 at 11:09):but state.ξₖ is a Float64. That's exactly the point that I don't get
Max Köhler (Jun 14 2021 at 11:10):Sebastian Pfitzner said:
ah, yeah, you'd want to annotate the rhs, not the lhs
so e.g. volumefraction(...)::Float64?
Fredrik Ekre (Jun 14 2021 at 11:13):Max Köhler said:
but
state.ξₖis a Float64. That's exactly the point that I don't get
Just to make sure, that field is strongly typed in the struct definition?
Max Köhler (Jun 14 2021 at 11:13):yes its
Base.@kwdef struct RelaxedDamageState{dim,T,M,N} <: MaterialState
damage::Damage{dim,T,M} = Damage()
damage⁺::Damage{dim,T,M} = Damage()
damage¯::Damage{dim,T,M} = Damage()
P::Tensor{2,dim,T,M} = Tensor{2,1,Float64,1}((0.0,))
σ::SymmetricTensor{2,dim,T,N} = SymmetricTensor{2,1,Float64,1}((0.0,))
ξₖ::Float64 = 0.0
dₖ::Float64 = 0.0
W_min::Float64 = 0.0
W::Float64 = 0.0
convex::Bool = true
relaxed::Bool = false
end
This is now a type stable version, but I still don't understand why. It's pretty much your suggestion @Andrey Oskin
function constitutive_driver(F, material::RelaxedDamage, state::RelaxedDamageState)
state.relaxed ? convex = false : convex = true
W_min = state.W_min
W = state.W
F⁺, F¯ = state.damage⁺.Fₖ, state.damage¯.Fₖ
#local ξ::Float64
#local d::Float64
ξ_old::Float64 = state.ξₖ
d_old::Float64 = state.dₖ
if convex
W_min, F⁺, F¯ = convexify(F, material, state)
#ξ_new = (F[1] - F¯[1]) / (F⁺[1] - F¯[1])::Float64
#d_new = (F⁺[1] - F¯[1]) / F[1]::Float64
W = W_energy(F, material, state.damage)
convex = W < W_min || isapprox(W_min,W)
end
if !(convex)
convex = false
relaxed = true
ξ_new = volumefraction(state, F, F⁺,F¯)::Float64
if (0.0 < ξ_new < 1.0) && (F[1] >= F⁺[1] || F[1] <= F¯[1])
relaxed = false
else
relaxed = true
end
d_new = (ξ_new != 1 ? (F⁺[1] - F¯[1]) / F[1] : 0.0)::Float64
#P, 𝔸, ψ⁺, ψ¯ = homogenized_derivatives(ξ, d, F, material, state) FIXME
𝔸, P = Tensors.hessian(
y -> W_homogenized(ξ_new, d_new, y, material, state),
F,
:all)
ψ⁺ = Ψ(tdot(F⁺), material.base_material)
ψ¯ = Ψ(tdot(F¯), material.base_material)
σ = symmetric(inv(det(F)) * P ⋅ F')
W = W_homogenized(ξ_new, d_new, F, material, state)
unrelaxed = Damage(ψₖ = state.damage.ψₖ, βₖ = state.damage.βₖ, Fₖ=F)
phase⁺ = Damage(ψₖ = ψ⁺, βₖ = max(state.damage⁺.βₖ, ψ⁺), Fₖ=F⁺)
phase¯ = Damage(ψₖ = ψ¯, βₖ = max(state.damage¯.βₖ, ψ¯), Fₖ=F¯)
newstate = RelaxedDamageState(unrelaxed,phase⁺,phase¯,P,σ, ξ_new,d_new, W_min, W, convex, relaxed)
else
#𝔸, P, ψ = strain_energy_derivatives(F, material, state)
𝔸, P = Tensors.hessian(y -> W_energy(y, material, state.damage), F, :all)
ψ = Ψ(tdot(F), material.base_material)
convex = true
relaxed = false
σ = symmetric(inv(det(F)) * state.P ⋅ F')
#ξ = copy(state.ξₖ)::Float64
#d = copy(state.dₖ)::Float64
unrelaxed = Damage(ψₖ = ψ, βₖ = max(state.damage.βₖ, ψ), Fₖ=F)
phase⁺ = Damage(ψₖ = ψ, βₖ = max(state.damage.βₖ, ψ))
phase¯ = Damage(ψₖ = ψ, βₖ = max(state.damage.βₖ, ψ))
newstate = RelaxedDamageState(unrelaxed,phase⁺,phase¯,P,σ, ξ_old,d_old, W_min, W, convex, relaxed)
end
return P, 𝔸, newstate
end
Hmm, it shouldn't work, since if convex is true, newstate is not defined.
Max Köhler (Jun 14 2021 at 11:25):Really? Both branches in the if-else statement define newstate
Max Köhler (Jun 14 2021 at 11:33):The type annotations are actually not needed in the above example. Unfortunately, I can't get a version working, which uses the local statement with type annotations. If I try to do it, I directly get a Core.Box at initialization.
Kwaku Oskin (Jun 14 2021 at 11:47):Ah, I see, yes, you are right. Yes, everything is defined properly in all branches.
Max Köhler (Jun 14 2021 at 12:23):So I achieved now somewhat what I wanted, i.e. a type stable version but I still don't understand why this fixes the Core.Boxproblem. I have some intuition what is meant by that, but in this particular case I cannot see the logic behind it. Other than that and quite a related question for a different function:
Is there a type stable version of findfirst or should I go for a small custom function?
Fredrik Ekre (Jun 14 2021 at 12:24):Just make sure to properly handle the nothing case and there shouldn't be a problem.
Max Köhler (Jun 14 2021 at 12:25):As far as I understood, GPUCompiler will complain about any type instability, but maybe the Union with nothing is an exception
Fredrik Ekre (Jun 14 2021 at 12:26):I mean, if you do something like
idx = findfirst(...)
idx === nothing && error(...)
# something with idx
isn't that okay?
Max Köhler (Jun 14 2021 at 12:31):ah ok nvm, findfirst seems to have it's own dispatch for gpuarrays, so should be okay I guess
Last updated: Nov 06 2024 at 04:40 UTC