Help with solving a DDE

Specific Domains Modelling & Simulations
1

Hi!

I am trying to use Julia’s DifferentialEquations.jl package to solve DDE. I was able to solve my problem for two objects using the code below

clearconsole()
@time using DifferentialEquations
const two_bubble =
  let σ=0.0725,  ρ=998,  μ=1e-3,  c=1481,  pvtinf=0,  pinf0=1.01e5,  k=7/5,
      dist=10e-6, τ=dist/c,  L=[Inf dist;dist Inf], x=[0 dist],
      ps=100e3,    f=1e6,    R=[1e-6, 1e-6],
      SP=0.01, cycle=5;

    global function two_bubble!(du, u, h, p, t)



    A1=(1+(1-3*k)*u[2]/c)*((pinf0-pvtinf)/ρ+2*σ/(ρ*R[1]))*(R[1]/u[1])^(3*k)-2*σ/(ρ*u[1])-4*μ*u[2]/(ρ*u[1])-(1+u[2]/c)*(pinf0-pvtinf+ps*sin(2*pi*f*t-(2*pi*f/c)*x[1]))/ρ-ps*u[1]*cos(2*pi*f*t-(2*pi*f/c)*x[1])*2*pi*f/(ρ*c);
    A2=(1+(1-3*k)*u[4]/c)*((pinf0-pvtinf)/ρ+2*σ/(ρ*R[2]))*(R[2]/u[3])^(3*k)-2*σ/(ρ*u[3])-4*μ*u[4]/(ρ*u[3])-(1+u[4]/c)*(pinf0-pvtinf+ps*sin(2*pi*f*t-(2*pi*f/c)*x[2]))/ρ-ps*u[3]*cos(2*pi*f*t-(2*pi*f/c)*x[2])*2*pi*f/(ρ*c);

    du[1]= u[2]
    du[2]= (A1-(1.5*(1-u[2]/(3*c)))*u[2]^2-(2*h(p, t - τ; idxs = 3)*h(p, t - τ; idxs = 4)^2-h(p, t - τ, Val{1}; idxs = 4)*h(p, t - τ; idxs = 3)^2)/L[1, 2])/((1-u[2]/c)*u[1]+4*μ/(ρ*c))
    du[3]= u[4]
    du[4]= (A2-(1.5*(1-u[4]/(3*c)))*u[4]^2-(2*h(p, t - τ; idxs = 1)*h(p, t - τ; idxs = 2)^2-h(p, t - τ, Val{1}; idxs = 2)*h(p, t - τ; idxs = 1)^2)/L[1, 2])/((1-u[4]/c)*u[3]+4*μ/(ρ*c));

      nothing
    end

    global function h_two_bubble(p, t; idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")
R=[1e-6, 1e-6]
      if idxs === nothing
        [R[1],0,R[2],0]
      elseif idxs == 1
      R[1]
      elseif idxs == 2
      0
      elseif idxs == 3
      R[2]
      elseif idxs == 4
      0
      else
        error("delay differential equation consists of two components")
      end
    end

    global function h_two_bubble(p, t, ::Type{Val{1}};
                                    idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      if idxs === nothing
        [0,0,0,0]
      elseif idxs == 1 || idxs == 3
        0
      elseif idxs == 2 || idxs==4
        0
      else
        error("delay differential equation consists of two components")
      end
    end


    DDEProblem(two_bubble!, h_two_bubble, (0.0, cycle/f);
               constant_lags = [τ], neutral = true)

end



f=1e6;
SP=0.01

@time u=  solve(two_bubble,Tsit5(),saveat=SP/f,span=SP/f)
r=convert(Array,u)
t=u.t

I understand that the written code is not optimized at all but at this point I am only trying to get it to work and work on optimization after.

clearconsole()


@time using DifferentialEquations
Sp=11; Sf=2.5e-6
SR0=0.01; kai=Sp/2; ks=Sf/(12*pi); RBB=1.02
      σ=0.0725;  ρ=998;  μ=1e-3;  c=1481;  pvtinf=0;  pinf0=1.01e5;  k=7/5;
      ps=1e3; f=1e6; R0=1.7e-6; Nb=2; d=100e-6; Mn=20e-6;
      SP=0.01; cycle=80;



R=zeros(1,Nb)
for i=1:Nb
  R[i]=R0
end
rbreak=R.*RBB; rb=R./sqrt(SR0/kai +1);


X=zeros(1,Nb)
Y=zeros(1,Nb)
Z=zeros(1,Nb)
L=zeros(Nb,Nb)
X[1]=d*rand(1)[1]
Y[1]=d*rand(1)[1]
Z[1]=d*rand(1)[1]
L[1,1]=Inf

k=1
while minimum(L)<Mn
global  i=Int64(2)
println("Attempt $k
")
while i<=Nb
X[i]=d*rand(1)[1]
Y[i]=d*rand(1)[1]
Z[i]=d*rand(1)[1]
for j=1:i
  if i==j
  L[i,j]=Inf
  else
  L[i,j]=sqrt(((X[i]-X[j])^2)+((Y[i]-Y[j])^2)+((Z[i]-Z[j])^2));
  L[j,i]=L[i,j];
  end
end
      global i=i+1
end
global k=k+1
end




τ=L/c;



p=[σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,rbreak,rb,X,τ]
################################################################################################################################################################################################################################################
function N_buckle!(du, u, h, p, t)  #with secondary delays
  σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,rbreak,rb,X,τ=p

  TT=zeros(1,Nb)
  P=zeros(1,Nb)
  SR=zeros(1,Nb)
  Psc=zeros(1,Nb)
  for i=1:Nb
   if 2*pi*f*t-(2*pi*f/c)*X[i]<0
        TT[i]=0
    else
        TT[i]=1
    end
    if u[2*i-1]>=rb[i] && u[2*i-1]<rbreak[i]
        SR[i]=kai*((u[2*i-1]/rb[i])^2-1)
    elseif u[2*i-1]>=rbreak[i]
        SR[i]=σ
    end

    P[i]=((pinf0+2*SR0)*(R[i]/u[2*i-1])^(3*k))*(1-3*k*u[2*i]/c)-(4*μ*u[2*i]/u[1])-(2*SR[i]/u[2*i-1])-(4*ks*u[2*i]/(u[2*i-1]^2))-pinf0-ps*TT[i]*sin(2*pi*f*t-(2*pi*f/c)*X[i]);


    Psc[i]=0
    for j=1:Nb
         if j~=i
         Pcs[i]=Pcs[i]+(2*h(p, t - τ[i,j]; idxs = 2*j-1)*h(p, t - τ[i,j]; idxs = 2*j)^2-h(p, t - τ[i,j], Val{1}; idxs = 2*j)*h(p, t - τ[i,j]; idxs = 2*j-1)^2)/(L[i,j])
         end
    end

du[2*i-1]=u[2*i]
du[2*i]=(P[2*i-1]/(ρ*u[2*i-1]))-(1.5*(u[2*i]^2)/u[2*i-1])-(Psc[i])/u[2*i-1]
end

  nothing
end
##################################################################################################
##################################################################################################
##################################################################################################
function h_N_buckle(p, t; idxs::Union{Nothing,Int} = nothing)
t ≤ 0 || error("history function is only implemented for t ≤ 0")
σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,rbreak,rb,X,τ=p

  u0=zeros(1,2*Nb)
  for i=1:2
  u0[2*i-1]=R[i]
  u0[2*i]=0
  end


if idxs==nothing
  u0
else
for i=1:2*Nb
 if idxs==2*i-1
  R[i]
elseif idxs==2*i
  0
 end
 end
end
end
##################################################################################################
##################################################################################################
##################################################################################################
function h_N_buckle(p, t, ::Type{Val{1}};
                                idxs::Union{Nothing,Int} = nothing)
t ≤ 0 || error("history function is only implemented for t ≤ 0")
σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,Rbreak,Rb,X,τ=p

  if idxs === nothing
    zeros(1,2*Nb)
     else
     for i=1:Nb
     0
     end
  end
end
##################################################################################################
##################################################################################################
##################################################################################################
MKMWD=DDEProblem(N_buckle!, h_N_buckle, (0.0, cycle/f),p;
               constant_lags = [τ], neutral = true)

SP=0.01;
tol=1e-4
alg=Tsit5()

u2=solve(MKMWD,alg,saveat=SP/f,reltol=tol,abstol=tol)

When I am trying to make my code (MultiDelayedBuckle.jl) scalable to solve for arbitrary number of objects I am encountering an error I am not able to understand

LoadError: MethodError: no method matching isless(::Matrix{Float64}, ::Float64)
Closest candidates are:
  isless(!Matched::Union{StatsBase.PValue, StatsBase.TestStat}, ::AbstractFloat) at C:\Users\hhagh\.julia\packages\StatsBase\Lc3YW\src\statmodels.jl:504
  isless(!Matched::Union{StatsBase.PValue, StatsBase.TestStat}, ::Real) at C:\Users\hhagh\.julia\packages\StatsBase\Lc3YW\src\statmodels.jl:495
  isless(!Matched::Discontinuity, ::Number) at C:\Users\hhagh\.julia\packages\DelayDiffEq\otmBn\src\discontinuity_type.jl:21
  ...
in expression starting at D:\Julia\Working\MultiDelayedBuckle.jl:155
<(x::Matrix{Float64}, y::Float64) at operators.jl:279
initialize_tstops_d_discontinuities_propagated(#unused#::Type{Float64}, tstops::Tuple{}, d_discontinuities::Tuple{}, tspan::Tuple{Float64, Float64}, order_discontinuity_t0::Int64, alg_maximum_order::Int64, constant_lags::Vector{Matrix{Float64}}, neutral::Bool) at solve.jl:438
__init(prob::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, alg::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}; saveat::Float64, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Float64, reltol::Float64, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, qoldinit::Rational{Int64}, fullnormalize::Bool, failfactor::Int64, beta1::Nothing, beta2::Nothing, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, discontinuity_interp_points::Int64, discontinuity_abstol::Float64, discontinuity_reltol::Int64, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:default_set,), Tuple{Bool}}}) at solve.jl:187
(::SciMLBase.var"#__init##kw")(::NamedTuple{(:default_set, :saveat, :reltol, :abstol), Tuple{Bool, Float64, Float64, Float64}}, ::typeof(SciMLBase.__init), prob::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, alg::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}) at solve.jl:67
__init at solve.jl:67 [inlined]
__solve(::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, ::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}; kwargs::Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:default_set, :saveat, :reltol, :abstol), Tuple{Bool, Float64, Float64, Float64}}}) at solve.jl:4
(::SciMLBase.var"#__solve##kw")(::NamedTuple{(:default_set, :saveat, :reltol, :abstol), Tuple{Bool, Float64, Float64, Float64}}, ::typeof(SciMLBase.__solve), ::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, ::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}) at solve.jl:4
__solve(::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, ::Tsit5; default_set::Bool, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :reltol, :abstol), Tuple{Float64, Float64, Float64}}}) at default_solve.jl:7
__solve at default_solve.jl:3 [inlined]
#solve_call#56 at solve.jl:61 [inlined]
solve_call at solve.jl:48 [inlined]
#solve_up#58 at solve.jl:82 [inlined]
solve_up at solve.jl:75 [inlined]
#solve#57 at solve.jl:70 [inlined]
(::CommonSolve.var"#solve##kw")(::NamedTuple{(:saveat, :reltol, :abstol), Tuple{Float64, Float64, Float64}}, ::typeof(solve), prob::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, args::Tsit5) at solve.jl:68
top-level scope at MultiDelayedBuckle.jl:155
eval at boot.jl:360 [inlined]
include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String) at loading.jl:1094

Any help would be greatly appreciated!

2

The error is pointing to line 155 and you shared a 147 line code, so this cannot be the same code you’re actually running.

My guess is what’s going on is that you needed .<.

P.S. Just do zeros(Nb) and simplify your code.

1 Like
4

You are right, here is the correct error message corresponding to correct line numbers

LoadError: MethodError: no method matching isless(::Matrix{Float64}, ::Float64)
Closest candidates are:
  isless(!Matched::Union{StatsBase.PValue, StatsBase.TestStat}, ::AbstractFloat) at C:\Users\hhagh\.julia\packages\StatsBase\Lc3YW\src\statmodels.jl:504
  isless(!Matched::Union{StatsBase.PValue, StatsBase.TestStat}, ::Real) at C:\Users\hhagh\.julia\packages\StatsBase\Lc3YW\src\statmodels.jl:495
  isless(!Matched::Discontinuity, ::Number) at C:\Users\hhagh\.julia\packages\DelayDiffEq\otmBn\src\discontinuity_type.jl:21
  ...
in expression starting at D:\Julia\Working\MultiDelayedBuckle.jl:147
<(x::Matrix{Float64}, y::Float64) at operators.jl:279
initialize_tstops_d_discontinuities_propagated(#unused#::Type{Float64}, tstops::Tuple{}, d_discontinuities::Tuple{}, tspan::Tuple{Float64, Float64}, order_discontinuity_t0::Int64, alg_maximum_order::Int64, constant_lags::Vector{Matrix{Float64}}, neutral::Bool) at solve.jl:438
__init(prob::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, alg::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}; saveat::Float64, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Float64, reltol::Float64, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, qoldinit::Rational{Int64}, fullnormalize::Bool, failfactor::Int64, beta1::Nothing, beta2::Nothing, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, discontinuity_interp_points::Int64, discontinuity_abstol::Float64, discontinuity_reltol::Int64, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:default_set,), Tuple{Bool}}}) at solve.jl:187
(::SciMLBase.var"#__init##kw")(::NamedTuple{(:default_set, :saveat, :reltol, :abstol), Tuple{Bool, Float64, Float64, Float64}}, ::typeof(SciMLBase.__init), prob::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, alg::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}) at solve.jl:67
__init at solve.jl:67 [inlined]
__solve(::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, ::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}; kwargs::Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:default_set, :saveat, :reltol, :abstol), Tuple{Bool, Float64, Float64, Float64}}}) at solve.jl:4
(::SciMLBase.var"#__solve##kw")(::NamedTuple{(:default_set, :saveat, :reltol, :abstol), Tuple{Bool, Float64, Float64, Float64}}, ::typeof(SciMLBase.__solve), ::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, ::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}) at solve.jl:4
__solve(::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, ::Tsit5; default_set::Bool, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :reltol, :abstol), Tuple{Float64, Float64, Float64}}}) at default_solve.jl:7
__solve at default_solve.jl:3 [inlined]
#solve_call#56 at solve.jl:61 [inlined]
solve_call at solve.jl:48 [inlined]
#solve_up#58 at solve.jl:82 [inlined]
solve_up at solve.jl:75 [inlined]
#solve#57 at solve.jl:70 [inlined]
(::CommonSolve.var"#solve##kw")(::NamedTuple{(:saveat, :reltol, :abstol), Tuple{Float64, Float64, Float64}}, ::typeof(solve), prob::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, args::Tsit5) at solve.jl:68
top-level scope at MultiDelayedBuckle.jl:147
eval at boot.jl:360 [inlined]
include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String) at loading.jl:1094
5

I tried this and it doesn’t work and now it returns the history function is only implemented for t ≤ 0 that is defined in the code. I also tried .≤ and it didn’t help either and resulted in the original error.

6

@ChrisRackauckas Do you have any other ideas that might be of help?

7

What’s the simplified code to copy paste and test?

8

This

clearconsole()


@time using DifferentialEquations
Sp=11; Sf=2.5e-6
SR0=0.01; kai=Sp/2; ks=Sf/(12*pi); RBB=1.02
      σ=0.0725;  ρ=998;  μ=1e-3;  c=1481;  pvtinf=0;  pinf0=1.01e5;  k=7/5;
      ps=1e3; f=1e6; R0=1.7e-6; Nb=2; d=100e-6; Mn=20e-6;
      SP=0.01; cycle=80;



R=zeros(1,Nb)
for i=1:Nb
  R[i]=R0
end
rbreak=R.*RBB; rb=R./sqrt(SR0/kai +1);


X=zeros(1,Nb)
Y=zeros(1,Nb)
Z=zeros(1,Nb)
L=zeros(Nb,Nb)
X[1]=d*rand(1)[1]
Y[1]=d*rand(1)[1]
Z[1]=d*rand(1)[1]
L[1,1]=Inf

k=1
while minimum(L)<Mn
global  i=Int64(2)
println("Attempt $k
")
while i<=Nb
X[i]=d*rand(1)[1]
Y[i]=d*rand(1)[1]
Z[i]=d*rand(1)[1]
for j=1:i
  if i==j
  L[i,j]=Inf
  else
  L[i,j]=sqrt(((X[i]-X[j])^2)+((Y[i]-Y[j])^2)+((Z[i]-Z[j])^2));
  L[j,i]=L[i,j];
  end
end
      global i=i+1
end
global k=k+1
end




τ=L/c;



p=[σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,rbreak,rb,X,τ]
################################################################################################################################################################################################################################################
function N_buckle!(du, u, h, p, t)  #with secondary delays
  σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,rbreak,rb,X,τ=p

  TT=zeros(1,Nb)
  P=zeros(1,Nb)
  SR=zeros(1,Nb)
  Psc=zeros(1,Nb)
  for i=1:Nb
   if 2*pi*f*t-(2*pi*f/c)*X[i]<0
        TT[i]=0
    else
        TT[i]=1
    end
    if u[2*i-1]>=rb[i] && u[2*i-1]<rbreak[i]
        SR[i]=kai*((u[2*i-1]/rb[i])^2-1)
    elseif u[2*i-1]>=rbreak[i]
        SR[i]=σ
    end

    P[i]=((pinf0+2*SR0)*(R[i]/u[2*i-1])^(3*k))*(1-3*k*u[2*i]/c)-(4*μ*u[2*i]/u[1])-(2*SR[i]/u[2*i-1])-(4*ks*u[2*i]/(u[2*i-1]^2))-pinf0-ps*TT[i]*sin(2*pi*f*t-(2*pi*f/c)*X[i]);


    Psc[i]=0
    for j=1:Nb
         if j~=i
         Pcs[i]=Pcs[i]+(2*h(p, t - τ[i,j]; idxs = 2*j-1)*h(p, t - τ[i,j]; idxs = 2*j)^2-h(p, t - τ[i,j], Val{1}; idxs = 2*j)*h(p, t - τ[i,j]; idxs = 2*j-1)^2)/(L[i,j])
         end
    end

du[2*i-1]=u[2*i]
du[2*i]=(P[2*i-1]/(ρ*u[2*i-1]))-(1.5*(u[2*i]^2)/u[2*i-1])-(Psc[i])/u[2*i-1]
end

  nothing
end
##################################################################################################
##################################################################################################
##################################################################################################
function h_N_buckle(p, t; idxs::Union{Nothing,Int} = nothing)
t ≤ 0 || error("history function is only implemented for t ≤ 0")
σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,rbreak,rb,X,τ=p

  u0=zeros(1,2*Nb)
  for i=1:2
  u0[2*i-1]=R[i]
  u0[2*i]=0
  end


if idxs==nothing
  u0
else
for i=1:2*Nb
 if idxs==2*i-1
  R[i]
elseif idxs==2*i
  0
 end
 end
end
end
##################################################################################################
##################################################################################################
##################################################################################################
function h_N_buckle(p, t, ::Type{Val{1}};
                                idxs::Union{Nothing,Int} = nothing)
t ≤ 0 || error("history function is only implemented for t ≤ 0")
σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,Rbreak,Rb,X,τ=p

  if idxs === nothing
    zeros(1,2*Nb)
     else
     for i=1:Nb
     0
     end
  end
end
##################################################################################################
##################################################################################################
##################################################################################################
MKMWD=DDEProblem(N_buckle!, h_N_buckle, (0.0, cycle/f),p;
               constant_lags = [τ], neutral = true)

SP=0.01;
tol=1e-4
alg=Tsit5()

u2=solve(MKMWD,alg,saveat=SP/f,reltol=tol,abstol=tol)
9

Can you simplify it a bit first? zeros(1,n) should be zeros(n). Etc. The tabbing is all off. That’s just really disrespectful of my time and I don’t have the time to do all of the cleanup before looking at the real issue. Does there really need to be all of those parameters in order to get the issue? You can’t delete a single parameter to recreate the issue? I don’t believe that.

2 Likes
10

oh I see, Yes that’s fair.

I managed to recreate the issue for the two object case and learned why it happens. just not sure how to address it
same error happens when I am addressing my delays in form of τ[1,2], τ[2,1] (see lines 36 and 38 where I am use them and line 5 where I am define them)

clearconsole()
@time using DifferentialEquations

      σ=0.0725;  ρ=998;  μ=1e-3;  c=1481;  pvtinf=0;  pinf0=1.01e5;  k=7/5;
      dist=10e-6;  L=[Inf dist;dist Inf]; τ=L/c; x=[0 dist];
      ps=1e3; f=5e6; R=[2e-6, 2e-6]; Nb=2;
      SR0=0.02; kai=0.9; ks=9e-9; RBB=1.01; rbreak=R.*RBB; rb=R./sqrt(SR0/kai +1);
      cycle=80;
      SP=0.01; tol=1e-5; alg=Tsit5()

################################################################################################################################################################################################################################################
    function two_buckle!(du, u, h, p, t) 
f=p[1];
      TT=zeros(Nb)
      P=zeros(Nb)

for i=1:Nb
      if 2*pi*f*t-(2*pi*f/c)*x[i]<0
            TT[i]=0
        else
            TT[i]=1
      end

        SR=zeros(1,2)
      if u[2*i-1]>=rb[i] && u[2*i-1]<rbreak[i]
            SR[i]=kai*((u[2*i-1]/rb[i])^2-1)
        elseif u[2*i-1]>=rbreak[i]
            SR[i]=σ
      end

        P[i]=((pinf0+2*SR0)*(R[i]/u[2*i-1])^(3*k))*(1-3*k*u[2*i]/c)-(4*μ*u[2*i]/u[1])-(2*SR[i]/u[2*i-1])-(4*ks*u[2*i]/(u[2*i-1]^2))-pinf0-ps*TT[i]*sin(2*pi*f*t-(2*pi*f/c)*x[i]);

end
   
    du[1]= u[2]
    du[2]= (P[1]/(ρ*u[1]))-(1.5*(u[2]^2)/u[1])-(2*h(p, t - τ[1,2]; idxs = 3)*h(p, t - τ[1,2]; idxs = 4)^2-h(p, t - τ[1,2], Val{1}; idxs = 4)*h(p, t - τ[1,2]; idxs = 3)^2)/(u[1]*L[1, 2])
    du[3]= u[4]
    du[4]= (P[2]/(ρ*u[3]))-(1.5*(u[4]^2)/u[3])-(2*h(p, t - τ[2,1]; idxs = 1)*h(p, t - τ[2,1]; idxs = 2)^2-h(p, t - τ[2,1], Val{1}; idxs = 2)*h(p, t - τ[2,1]; idxs = 1)^2)/(u[3]*L[2, 1])

      nothing
    end
########################################################################################################################################################################################################################################################################################
    u0=zeros(2*Nb)
    for i=1:2
    u0[2*i-1]=R[i]
    u0[2*i]=0
    end
    function h_two_buckle(p, t; idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      if idxs === nothing
      u0
      elseif mod(idxs,2)==1
      R[Int64((idxs+1)/2)]
      elseif mod(idxs,2)==0
      0
      end
    end
########################################################################################################################################################################################################################################################################################
   function h_two_buckle(p, t, ::Type{Val{1}};
                                    idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

       if idxs === nothing
      zeros(Nb)
    elseif mod(idxs,2)==1
    0
    elseif mod(idxs,2)==0
    0
      end
    end
########################################################################################################################################################################################################################################################################################
p=[f]
MKMWD=DDEProblem(two_buckle!, h_two_buckle, (0.0, cycle/f),p;
               constant_lags = [τ], neutral = true)


u=  DifferentialEquations.solve(MKMWD,alg,saveat=SP/f,reltol=tol,abstol=tol)

r=convert(Array,u)
t=u.t

but when I use the delays like below (see lines 5, 36 and 38) it runs with no issues

clearconsole()
@time using DifferentialEquations

      σ=0.0725;  ρ=998;  μ=1e-3;  c=1481;  pvtinf=0;  pinf0=1.01e5;  k=7/5;
      dist=10e-6;  L=[Inf dist;dist Inf]; τ=dist/c; x=[0 dist];
      ps=1e3; f=5e6; R=[2e-6, 2e-6]; Nb=2;
      SR0=0.02; kai=0.9; ks=9e-9; RBB=1.01; rbreak=R.*RBB; rb=R./sqrt(SR0/kai +1);
      cycle=80;
      SP=0.01; tol=1e-5; alg=Tsit5()

################################################################################################################################################################################################################################################
    function two_buckle!(du, u, h, p, t) 
f=p[1];
      TT=zeros(Nb)
      P=zeros(Nb)

for i=1:Nb
      if 2*pi*f*t-(2*pi*f/c)*x[i]<0
            TT[i]=0
        else
            TT[i]=1
      end

        SR=zeros(1,2)
      if u[2*i-1]>=rb[i] && u[2*i-1]<rbreak[i]
            SR[i]=kai*((u[2*i-1]/rb[i])^2-1)
        elseif u[2*i-1]>=rbreak[i]
            SR[i]=σ
      end

        P[i]=((pinf0+2*SR0)*(R[i]/u[2*i-1])^(3*k))*(1-3*k*u[2*i]/c)-(4*μ*u[2*i]/u[1])-(2*SR[i]/u[2*i-1])-(4*ks*u[2*i]/(u[2*i-1]^2))-pinf0-ps*TT[i]*sin(2*pi*f*t-(2*pi*f/c)*x[i]);

end
   
    du[1]= u[2]
    du[2]= (P[1]/(ρ*u[1]))-(1.5*(u[2]^2)/u[1])-(2*h(p, t - τ; idxs = 3)*h(p, t - τ; idxs = 4)^2-h(p, t - τ, Val{1}; idxs = 4)*h(p, t - τ; idxs = 3)^2)/(u[1]*L[1, 2])
    du[3]= u[4]
    du[4]= (P[2]/(ρ*u[3]))-(1.5*(u[4]^2)/u[3])-(2*h(p, t - τ; idxs = 1)*h(p, t - τ; idxs = 2)^2-h(p, t - τ, Val{1}; idxs = 2)*h(p, t - τ; idxs = 1)^2)/(u[3]*L[2, 1])

      nothing
    end
########################################################################################################################################################################################################################################################################################
    u0=zeros(2*Nb)
    for i=1:2
    u0[2*i-1]=R[i]
    u0[2*i]=0
    end
    function h_two_buckle(p, t; idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      if idxs === nothing
      u0
      elseif mod(idxs,2)==1
      R[Int64((idxs+1)/2)]
      elseif mod(idxs,2)==0
      0
      end
    end
########################################################################################################################################################################################################################################################################################
   function h_two_buckle(p, t, ::Type{Val{1}};
                                    idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

       if idxs === nothing
      zeros(Nb)
    elseif mod(idxs,2)==1
    0
    elseif mod(idxs,2)==0
    0
      end
    end
########################################################################################################################################################################################################################################################################################
p=[f]
MKMWD=DDEProblem(two_buckle!, h_two_buckle, (0.0, cycle/f),p;
               constant_lags = [τ], neutral = true)


u=  DifferentialEquations.solve(MKMWD,alg,saveat=SP/f,reltol=tol,abstol=tol)

r=convert(Array,u)
t=u.t

This is possible for the case of two objects in my equations since τ[1,2]=τ[2,1] but I cannot generalize it for the case of N objects. I also need to add the τ[i,j]s are constant and do not change throughout my solution.

11

Sorry it took a bit to get back to this. The answer is rather simple:

MKMWD=DDEProblem(two_buckle!, h_two_buckle, (0.0, cycle/f),p;
               constant_lags = τ, neutral = true)

You had constant_lags = [τ]. constant_lags is supposed to be an array of constant lag times. Your τ is already an array:

c=1481
L=[Inf dist;dist Inf]; τ=L/c

so putting that array inside of another array confused it, and it was trying to act on τ as though it was a scalar. We could probably throw a better error for that, but hopefully that gets you working.

12

Thank you this solved my issue.