`jacobian_ad_forward` Not working with `IndicatorHennemannGassner`

[DanielDoehring]

Not sure if this is a known issue, but right now calling jacobian_ad_forward with a semidiscretization build from VolumeIntegralShockCapturingHG with shock indicator IndicatorHennemannGassner like in examples/structured_1d_dgsem/elixir_euler_sedov.jl gives the error

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{Nothing, Float64, 12})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:772
  (::Type{T})(::VectorizationBase.Double{T}) where T<:Union{Float16, Float32, Float64, VectorizationBase.Vec{<:Any, <:Union{Float16, Float32, Float64}}, VectorizationBase.VecUnroll{var"#s36", var"#s35", var"#s34", V} where {var"#s36", var"#s35", var"#s34"<:Union{Float16, Float32, Float64}, V<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8, SIMDTypes.Bit, VectorizationBase.AbstractSIMD{var"#s35", var"#s34"}}}} at ~/.julia/packages/VectorizationBase/MOxkH/src/special/double.jl:100
  ...
Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{Nothing, Float64, 12})
    @ Base ./number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{Nothing, Float64, 12}, i1::Int64)
    @ Base ./array.jl:966
  [3] macro expansion
    @ ~/.julia/dev/Trixi/src/solvers/dgsem_tree/indicators_1d.jl:51 [inlined]

...

i.e., what goes wrong is that the datatypes within

indicator[i] = indicator_hg.variable(u_local, equations)

disagree.

[sloede]

Hm. Can you you add

@show typeof(indictor[i])
@show typeof(u_local)
@show typeof(indicator_hg.variable)

and post the output? AFAICT, u_local should be correctly typed as a Dual (it comes from get_node_vars) andindicator_hg.variable should be density_pressure, which also looks innocent

However, indicator[i] is of type real(basis) and I suspect that this might be a Float64.

[ranocha]

This is already tracked in https://github.com/trixi-framework/Trixi.jl/issues/462. You can find a reference to some description of the issue there. Of course, it would be great to fix this. Maybe https://github.com/SciML/PreallocationTools.jl could be worth a try?

[DanielDoehring]

This is already tracked in #462. You can find a reference to some description of the issue there. Of course, it would be great to fix this. Maybe https://github.com/SciML/PreallocationTools.jl could be worth a try?

Ah I see, I couldn't find the issue with the search function - I will close this for now to avoid duplicate issues.

[ranocha]

No, let's keep this one open to track this specific request - would be great to get someone working on this :slightly_smiling_face:

[DanielDoehring]

typeof(indicator[i]) = Float64
typeof(u_local) = SVector{3, ForwardDiff.Dual{Nothing, Float64, 12}}
typeof(indicator_hg.variable) = typeof(density_pressure)
typeof(indicator_hg.variable(u_local, equations)) = ForwardDiff.Dual{Nothing, Float64, 12}

[DanielDoehring]

I found the time to look into this again.

Using indicator[i] = indicator_hg.variable(u_local, equations).value instead of indicator[i] = indicator_hg.variable(u_local, equations) seems to solve the issue. Truncating the DualNumber to its value stops of course the further propagation of derivatives but this seems to work out (at least now) as there are no further errors, i.e., the indicator array seems not to be expected to carry derivative information.

Computing the spectrum of the semidiscretization of the tree_1d_dgsem/elixir_euler_blast_wave.jl example and comparing to the spectrum of a basic DGSEM solver DGSEM(basis, surface_flux) reveals that the spectra at least differ.

[ranocha]

By doing this, you basically compute the spectrum of the blended RHS for fixed indicator values. If this is what you need right now - great! A "real" solution should allow us to differentiate through the indicator completely.

[DanielDoehring]

For the moment that might be actually enough, yes.

Nevertheless, I kept tinkering around for a bit longer - for the sake of documentation I poste that here. So the reason why the types in https://github.com/trixi-framework/Trixi.jl/blob/3ea50d1d108d9c3b4d5ce463a40a20bfa7937633/src/solvers/dgsem_tree/indicators_1d.jl#L51 inherently disagree is that u_local gets promoted to a type with ForwardDiff.Dual as a base in the call of jacobian_ad_forward: https://github.com/trixi-framework/Trixi.jl/blob/7dd423222aa800c635bbdbdcb3c3211736e3d601/src/semidiscretization/semidiscretization.jl#L239

Although the solver gets remade, https://github.com/trixi-framework/Trixi.jl/blob/7dd423222aa800c635bbdbdcb3c3211736e3d601/src/semidiscretization/semidiscretization_hyperbolic.jl#L78

this seems not to propagate through to the indicator, which still bears the datatype assigned at creation:

https://github.com/trixi-framework/Trixi.jl/blob/7dd423222aa800c635bbdbdcb3c3211736e3d601/src/solvers/dgsem_tree/indicators_1d.jl#L15

I tried to manually promote the type to ForwardDiff types by adding indicator = convert(Vector{typeof(indicator_hg.variable(get_node_vars(u, equations, dg, 1, element), equations))}, indicator)

after

https://github.com/trixi-framework/Trixi.jl/blob/7dd423222aa800c635bbdbdcb3c3211736e3d601/src/solvers/dgsem_tree/indicators_1d.jl#L45

and substitute

   Axes = axes(data_out, 1)
  data_out = convert(typeof(data_in), data_out)

  # @tullio threads=false data_out[i] = matrix[i, ii] * data_in[ii]
  #@turbo for i in axes(data_out, 1)
    for i in Axes

for

https://github.com/trixi-framework/Trixi.jl/blob/7dd423222aa800c635bbdbdcb3c3211736e3d601/src/solvers/dgsem/interpolation.jl#L129

When doing so, one can obtain a Jacobian - with is full of NaN. I have the feeling that this might be related to the initialization of the ForwardDiff type https://github.com/trixi-framework/Trixi.jl/blob/7dd423222aa800c635bbdbdcb3c3211736e3d601/src/semidiscretization/semidiscretization.jl#L230

which would probably needed to be extended by the parameters of the indicator for this to work out.