Merge pull request #346 from CliMA/ck/zeros

Use zero instead of similar to avoid undefined behavior

Author: charleskawczynski

src/nl_solvers/newtons_method.jl

@@ -149,7 +149,7 @@ Base.@kwdef struct ForwardDiffJVP{S <: ForwardDiffStepSize, T} <: JacobianFreeJV
     step_adjustment::T = 1
 end
 
-allocate_cache(::ForwardDiffJVP, x_prototype) = (; x2 = similar(x_prototype), f2 = similar(x_prototype))
+allocate_cache(::ForwardDiffJVP, x_prototype) = (; x2 = zero(x_prototype), f2 = zero(x_prototype))
 
 function jvp!(alg::ForwardDiffJVP, cache, jΔx, Δx, x, f!, f, post_implicit!)
     (; default_step, step_adjustment) = alg
@@ -562,8 +562,8 @@ function allocate_cache(alg::NewtonsMethod, x_prototype, j_prototype = nothing)
         convergence_checker_cache = isnothing(convergence_checker) ? nothing :
                                     allocate_cache(convergence_checker, x_prototype),
         Δx = similar(x_prototype),
-        f = similar(x_prototype),
-        j = isnothing(j_prototype) ? nothing : similar(j_prototype),
+        f = zero(x_prototype),
+        j = isnothing(j_prototype) ? nothing : zero(j_prototype),
     )
 end
 

src/solvers/imex_ssprk.jl

@@ -19,13 +19,13 @@ function init_cache(prob::DiffEqBase.AbstractODEProblem, alg::IMEXAlgorithm{SSP}
     s = length(b_exp)
     inds = ntuple(i -> i, s)
     inds_T_imp = filter(i -> !all(iszero, a_imp[:, i]) || !iszero(b_imp[i]), inds)
-    U = similar(u0)
-    U_exp = similar(u0)
-    T_lim = similar(u0)
-    T_exp = similar(u0)
-    U_lim = similar(u0)
-    T_imp = SparseContainer(map(i -> similar(u0), collect(1:length(inds_T_imp))), inds_T_imp)
-    temp = similar(u0)
+    U = zero(u0)
+    U_exp = zero(u0)
+    T_lim = zero(u0)
+    T_exp = zero(u0)
+    U_lim = zero(u0)
+    T_imp = SparseContainer(map(i -> zero(u0), collect(1:length(inds_T_imp))), inds_T_imp)
+    temp = zero(u0)
     â_exp = SparseCoeffs(vcat(a_exp.coeffs, b_exp.coeffs'))
     β = SparseCoeffs(diag(â_exp, -1))
     for i in 1:length(β)

src/solvers/mis.jl

@@ -63,12 +63,12 @@ function init_cache(
 
     tab = tableau(alg, eltype(prob.u0))
     N = length(tab.d)
-    ΔU = ntuple(n -> similar(prob.u0), N)
+    ΔU = ntuple(n -> zero(prob.u0), N)
     # at time i, contains
     #   ΔU[j] = U[j] - u    j < i
     #   ΔU[i] = U[i]
     #   ΔU[N] = offset vec
-    F = ntuple(n -> similar(prob.u0), N)
+    F = ntuple(n -> zero(prob.u0), N)
     return MultirateInfinitesimalStepCache(ΔU, F, tab)
 end
 

src/solvers/rosenbrock.jl

@@ -94,7 +94,7 @@ function init_cache(prob::DiffEqBase.AbstractODEProblem, alg::RosenbrockAlgorith
     fU_exp = zero(prob.u0)
     fU_lim = zero(prob.u0)
     ∂Y∂t = zero(prob.u0)
-    k = ntuple(n -> similar(prob.u0), Nstages)
+    k = ntuple(n -> zero(prob.u0), Nstages)
     if !isnothing(prob.f.T_imp!)
         W = prob.f.T_imp!.jac_prototype
     else

src/solvers/wickerskamarock.jl

@@ -26,8 +26,8 @@ struct WickerSkamarockRungeKuttaCache{T <: WickerSkamarockRungeKuttaTableau, A}
     F::A
 end
 function init_cache(prob::DiffEqBase.ODEProblem, alg::WickerSkamarockRungeKutta; kwargs...)
-    U = similar(prob.u0)
-    F = similar(prob.u0)
+    U = zero(prob.u0)
+    F = zero(prob.u0)
     return WickerSkamarockRungeKuttaCache(tableau(alg, eltype(F)), U, F)
 end
 

test/convergence_utils.jl

@@ -12,7 +12,7 @@ DirectSolver(args...) = DirectSolver()
 
 function (::DirectSolver)(x, A, b, matrix_updated; kwargs...)
     n = length(x)
-    M = mapslices(y -> mul!(similar(y), A, y), Matrix{eltype(x)}(I, n, n), dims = 1)
+    M = mapslices(y -> mul!(zero(y), A, y), Matrix{eltype(x)}(I, n, n), dims = 1)
     x .= M \ b
 end