Help with solving 2nd order PDE

David_F · September 20, 2022, 8:14am

Dear community,

This is my first attempt with Julia, so naturally I already face several problems, but I’m very excited to learn something new!

Basically, I’m trying to solve a 2nd order PDE, namely an extended Rayleigh-Plesset equation, which takes into account material strength by adding two stress terms, Pᵥ and Pᵧ, in addition to the standard Pₛ (external time-dependent pressure pulse) plus a (constant) porosity term. Both Pᵥ and Pᵧ depend on the bubble radius (denoted by a and/or its time derivative, da).

Can anybody take a look at the code I wrote below, and point to any immediate errors? For example, it seems that the parameter c=c_radius(a, b_radius(a)) given to the solver should be defined differently, because of “a is undefined” error. In addition, I’m not sure how to properly implement the if-else-if statements to adapt the stresses to the current value of α (which depends on the radius, a).

Thanks in advance!

using DifferentialEquations, Plots                                          
                                                                                                                                                                            
#Parameters                                                     
Y = 2.5e9                       #Yield stress (N/m2)                       
µ = 0.2e9                       #shear modulus (N/m2)                     
ρ = 1000.0                      #density (Kg.m-3)                 
η = 1e2                         #dynamic viscosity ((N/m2).s)                 
                                                                                                            
#Porosity
φ = 0.1

#Initial conditions                                                              
a₀ = 10e-6                      # bubble inner radius (m)
b₀ = a₀/φ^(1/3)                 # bubble outer radius (m)                           
da₀ = 1.0e-2                    #initial radial velocity (ms-1)            
tspan = (0.0, 100000.0e-6)      #(s) -- 1 acoustic cycle = 50 μs         

#Distention ratios
α₀ = b₀^3/(b₀^3 - a₀^3)
α₁ = (2*µ*α₀ + Y)/(2*µ+Y)
α₂ = 2*µ*α₀/(2*µ+Y)

# external drive
Pₛ = t->1e5*sin(2*pi*t)

# function definitions
function α(a, b)
   b^3/(b^3 - a^3) 
end

function b_radius(a)
   a/φ^(1/3) 
end

function c_radius(α)
   (2*µ/Y*a₀^(1/3)*(α₀-α)/(α₀-1))^(1/3)
end

function Pᵥ(a, c, da)
    if (α(a, b_radius(a)) <= α₀) && (α(a, b_radius(a)) >= α₁)
        0
    elseif (α(a, b_radius(a)) <= α₁) && (α(a, b_radius(a)) >= α₂)
        12*a*a*da*η*(-(1/3)*((1/c - 1/a)))
    elseif (α(a, b_radius(a)) <= α₂) && (α(a, b_radius(a)) >= 1)
        12*a*a*da*η*(-(1/3)*((1/b_radius(a)) - 1/a))        
    else
        0
    end
end

function Pᵧ(a)
    if (α(a, b_radius(a)) <= α₀) && (α(a, b_radius(a)) >= α₁)
        4*µ*(α₀-α(a,b_radius(a)))/(3*α(a,b_radius(a))*(α(a,b_radius(a))-1))
    elseif (α(a, b_radius(a)) <= α₁) && (α(a, b_radius(a)) >= α₂)
        (2*Y/3)*(1-(2*µ)/(Y*α(a,b_radius(a))) + log((2*µ*(α₀-α(a,b_radius(a))))/(Y*(α(a,b_radius(a))-1))))
    elseif (α(a, b_radius(a)) <= α₂) && (α(a, b_radius(a)) >= 1)
        (2*Y/3)*log(α(a,b_radius(a))/(α(a,b_radius(a))-1))
    else
        0
    end
end

function rpnnp(da, a, p, c, t)
    @. (((p(t)+Pᵧ(a)+Pᵥ(a, c, da))/ρ)-0.5*da*(4*da*(1-φ^(1/3))-da*(1-φ^(4/3))))/(a*(1-φ^(1/3)))
end


problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, Pₛ, c=c_radius(a, b_radius(a)))
solution = solve(problem, DPRKN6())

plot(solution, idxs=(2), linewidth=2, xaxis = "t ", yaxis = "R", framestyle = :box)

ChrisRackauckas · September 20, 2022, 8:32am

This is not valid syntax. You will have gotten a big warning that c is not a valid keyword argument. Heed that warning. The parameters should be one struct.

This is also not valid syntax. If you are writing for the non-mutating form, the function should be rpnnp(da, a, p, t). There is only one spot for parameters. Put your parameters into a single struct.

David_F · September 20, 2022, 10:07am

Thanks, Chris!
Since c depends on a (i.e the solution), I can’t define a struct outside the rpnnp function. So, as a workaround I just pass in the Pₛ variable, and let the c variable to be computed inside the rpnnp function.

It seems that the original problem is solved. However, there’s a new error, which I think is unrelated, namely:

MethodError: no method matching +(::var"#1#2", ::Float64)

I’m not sure which line it is complaining about, and what’s the cause for this.

The new code:


using DifferentialEquations, Plots                                          
                                                                                                                                                                            
#Parameters                                                     
Y = 2.5e9                       #Yield stress (N/m2)                       
µ = 0.2e9                       #shear modulus (N/m2)                     
ρ = 1000.0                      #density (Kg.m-3)                 
η = 1e2                         #dynamic viscosity ((N/m2).s)                 
                                                                                                            
#Porosity
φ = 0.1

#Initial conditions                                                              
a₀ = 10e-6                      # bubble inner radius (m)
b₀ = a₀/φ^(1/3)                 # bubble outer radius (m)                           
da₀ = 1.0e-2                    #initial radial velocity (ms-1)            
tspan = (0.0, 100000.0e-6)      #(s) -- 1 acoustic cycle = 50 μs         

#Distention ratios
α₀ = b₀^3/(b₀^3 - a₀^3)
α₁ = (2*µ*α₀ + Y)/(2*µ+Y)
α₂ = 2*µ*α₀/(2*µ+Y)

# external drive
Pₛ = t->1e5*sin(2*pi*t)

# function definitions
function α(a, b)
   b^3/(b^3 - a^3) 
end

function b_radius(a)
   a/φ^(1/3) 
end

function c_radius(α)
   (2*µ/Y*a₀^(1/3)*(α₀-α)/(α₀-1))^(1/3)
end

function Pᵥ(a, c, da)
    if (α(a, b_radius(a)) <= α₀) && (α(a, b_radius(a)) >= α₁)
        0
    elseif (α(a, b_radius(a)) <= α₁) && (α(a, b_radius(a)) >= α₂)
        12*a*a*da*η*(-(1/3)*((1/c - 1/a)))
    elseif (α(a, b_radius(a)) <= α₂) && (α(a, b_radius(a)) >= 1)
        12*a*a*da*η*(-(1/3)*((1/b_radius(a)) - 1/a))        
    else
        0
    end
end

function Pᵧ(a)
    if (α(a, b_radius(a)) <= α₀) && (α(a, b_radius(a)) >= α₁)
        4*µ*(α₀-α(a,b_radius(a)))/(3*α(a,b_radius(a))*(α(a,b_radius(a))-1))
    elseif (α(a, b_radius(a)) <= α₁) && (α(a, b_radius(a)) >= α₂)
        (2*Y/3)*(1-(2*µ)/(Y*α(a,b_radius(a))) + log((2*µ*(α₀-α(a,b_radius(a))))/(Y*(α(a,b_radius(a))-1))))
    elseif (α(a, b_radius(a)) <= α₂) && (α(a, b_radius(a)) >= 1)
        (2*Y/3)*log(α(a,b_radius(a))/(α(a,b_radius(a))-1))
    else
        0
    end
end

function rpnnp(da, a, Pₛ, t)
    c = c_radius(α(a, b_radius(a)))
    @. (((Pₛ+Pᵧ(a)+Pᵥ(a, c, da))/ρ)-0.5*da*(4*da*(1-φ^(1/3))-da*(1-φ^(4/3))))/(a*(1-φ^(1/3)))
end


problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, Pₛ)
solution = solve(problem, DPRKN6())

plot(solution, idxs=(2), linewidth=2, xaxis = "t ", yaxis = "R", framestyle = :box)

David_F · September 20, 2022, 11:34am

OK, I think I’m making progress…
I replaced the line

function rpnnp(da, a, Pₛ, t)

with

function rpnnp(da, a, p, t)

and p with p(t) in the line:

 @. (((p(t)+Pᵧ(a)+Pᵥ(a, c, da))/ρ)-0.5*da*(4*da*(1-φ^(1/3))-da*(1-φ^(4/3))))/(a*(1-φ^(1/3)))

and the error now seems of a better nature:

DomainError with -6.88868063402273e-18:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

I think the above error is easier to fix

Oscar_Smith · September 20, 2022, 12:03pm

use cbrt instead of x^(1/3). it will be faster and work better for negative inputs.

David_F · September 20, 2022, 12:12pm

Thanks, Oscar!
It seems that I get a solution, but only up to some point where the solution breaks due to dt being too small. So, I replaced DPRKN6() with Rodas4(), but then this error came up:

MethodError: no method matching vec(::Float64)
Closest candidates are:
  vec(::Distributions.Distribution{<:Distributions.ArrayLikeVariate}) at C:\Users\DAVIDFUR\.julia\packages\Distributions\39PV5\src\reshaped.jl:143
  vec(::ChainRulesCore.AbstractThunk) at C:\Users\DAVIDFUR\.julia\packages\ChainRulesCore\Z4Jry\src\tangent_types\thunks.jl:52
  vec(::StrideArraysCore.StrideArray) at C:\Users\DAVIDFUR\.julia\packages\StrideArraysCore\VQxXL\src\reshape.jl:11
  ...

Any idea what’s going on?

Oscar_Smith · September 20, 2022, 1:10pm

Can you post the stack trace? That has the information about where the error is.

David_F · September 20, 2022, 1:28pm

Ofc, here’s the stack trace

MethodError: no method matching vec(::Float64)
Closest candidates are:
vec(::Distributions.Distribution{<:Distributions.ArrayLikeVariate}) at C:\Users\DAVIDFUR.julia\packages\Distributions\39PV5\src\reshaped.jl:143
vec(::ChainRulesCore.AbstractThunk) at C:\Users\DAVIDFUR.julia\packages\ChainRulesCore\Z4Jry\src\tangent_types\thunks.jl:52
vec(::StrideArraysCore.StrideArray) at C:\Users\DAVIDFUR.julia\packages\StrideArraysCore\VQxXL\src\reshape.jl:11
…

Stacktrace:
[1] copyto!(dest::Vector{Float64}, A::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ RecursiveArrayTools C:\Users\DAVIDFUR.julia\packages\RecursiveArrayTools\IgoeN\src\array_partition.jl:171
[2] copyto_axcheck!(dest::Vector{Float64}, src::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ Base .\abstractarray.jl:1127
[3] Vector{Float64}(x::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ Base .\array.jl:626
[4] (Vector)(x::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ Core .\boot.jl:481
[5] zeromatrix(A::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ RecursiveArrayTools C:\Users\DAVIDFUR.julia\packages\RecursiveArrayTools\IgoeN\src\array_partition.jl:371
[6] build_J_W(alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, u::ArrayPartition{Float64, Tuple{Float64, Float64}}, uprev::ArrayPartition{Float64, Tuple{Float64, Float64}}, p::Function, t::Float64, dt::Float64, f::DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, #unused#::Type{Float64}, #unused#::Val{false})
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\derivative_utils.jl:878
[7] alg_cache(alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, u::ArrayPartition{Float64, Tuple{Float64, Float64}}, rate_prototype::ArrayPartition{Float64, Tuple{Float64, Float64}}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, uprev::ArrayPartition{Float64, Tuple{Float64, Float64}}, uprev2::ArrayPartition{Float64, Tuple{Float64, Float64}}, f::DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, t::Float64, dt::Float64, reltol::Float64, p::var"#53#54", calck::Bool, #unused#::Val{false})
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\caches\rosenbrock_caches.jl:506
[8] __init(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, 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::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Rational{Int64}, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, 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, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\solve.jl:312
[9] __init(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}) (repeats 5 times)
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\solve.jl:10
[10] __solve(::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, ::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\solve.jl:5
[11] __solve(::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, ::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing})
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\solve.jl:1
[12] solve_call(_prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; merge_callbacks::Bool, kwargshandle::KeywordArgError, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:470
[13] solve_call(_prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:440
[14] solve_up(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, sensealg::Nothing, u0::ArrayPartition{Float64, Tuple{Float64, Float64}}, p::Function, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:825
[15] solve_up(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, sensealg::Nothing, u0::ArrayPartition{Float64, Tuple{Float64, Float64}}, p::Function, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:798
[16] solve(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:795
[17] solve(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:786
[18] top-level scope
@ In[27]:72
[19] eval
@ .\boot.jl:368 [inlined]
[20] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1428

ChrisRackauckas · September 20, 2022, 10:32pm

Put your parameters in a tuple, and use the tuple.

params = (p, c)
problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, params)

function rpnnp(da, a, params, t)
     p, c = params
    @. (((p(t)+Pᵧ(a)+Pᵥ(a, c, da))/ρ)-0.5*da*(4*da*(1-φ^(1/3))-da*(1-φ^(4/3))))/(a*(1-φ^(1/3)))
end

David_F · September 29, 2022, 7:56am

Thanks, Chris!

After taking some time off from the digital world, I’m finally back.
I couldn’t use your suggestion because of an Undef Error - namely, c is a function of a, which is the solution, and so should be updated every iteration. If I define it as you typed it, I get an undefined error on a.
Actually, this leads me to my next question: how do I make sure the functions Pᵥ and Pᵧ update themselves according to the new value of a?

I need Pᵥ and Pᵧ (and c) to be updated at every iteration of the solution, but I’m not sure how to implement it.

Perhaps you or anybody else can assist?

Here’s the full, updated code:

using DifferentialEquations, Plots                                          
                                                                                                                                                                            
#Parameters                                                     
Y = 2.0e7                       #Yield stress (N/m2)                       
µ = 4.5e9                       #shear modulus (N/m2)                     
ρ = 1840.0                      #density (Kg.m-3)                 
η = 1.0e2                       #dynamic viscosity ((N/m2).s)                 
                                                                                                            
#Porosity
φ = 0.01

#Initial conditions                                                              
a₀ = 10e-6                      # bubble inner radius (m)
b₀ = a₀/cbrt(φ)                 # bubble outer radius (m)
da₀ = 1e-9                      # initial radial velocity (ms-1)            
tspan = (0.0, 5.0e-4)           # time (s)         

#Distention ratios
α₀ = b₀^3/(b₀^3 - a₀^3)
α₁ = (2*µ*α₀ + Y)/(2*µ+Y)
α₂ = 2*µ*α₀/(2*µ+Y)
print("α₀ = ",α₀,"  α₁=",α₁,"  α₂=",α₂)

Pₛ = t->1.0e16*t

# function definitions
function α(a, b)
   b^3/(b^3 - a^3)   
end

function b_radius(a)
   a/cbrt(φ)
end

function c_radius(α)
   cbrt((2*µ*(a₀^3)*(α₀-α))/(Y*(α₀-1)))
end

function Pᵥ(a, c, da, new_α)
    if (new_α <= α₀) && (new_α >= α₁)
        0
    elseif (new_α <= α₁) && (new_α >= α₂)
        12*a*a*da*η*(-(1/3)*((1/c^3 - 1/a^3)))
    elseif (new_α <= α₂) && (new_α >= 1)
        12*a*a*da*η*(-(1/3)*((1/b_radius(a)^3) - 1/a^3))        
    else
        0
    end
end

function Pᵧ(a, new_α)
    if (new_α <= α₀) && (new_α >= α₁)
        4*µ*(α₀-new_α)/(3*new_α*(new_α-1))
    elseif (new_α <= α₁) && (new_α >= α₂)
        (2*Y/3)*((1-(2*µ)/(Y*new_α) + log((2*µ*(α₀-new_α))/(Y*(new_α-1)))))
    elseif (new_α <= α₂) && (new_α >= 1)
        (2*Y/3)*log(new_α/(new_α-1))
    else
        0
    end
end

function rpnnp(da, a, p, t)
    new_α = α(a, b_radius(a))
    c = c_radius(α(a, b_radius(a)))
    @. (((-p(t)+Pᵧ(a, new_α)-Pᵥ(a, c, da, new_α))/ρ)-0.5*da*(4*da*(1-cbrt(φ))-da*(1-cbrt(φ)^4)))/(a*(1-cbrt(φ)))
end

function condition(y, t, integrator)
   #rise time of shock
   t >= 1.0e-17
end

function affect!(integrator)
   integrator.p[1] = 1.0e-9
end
cb = ContinuousCallback(condition, affect!)

problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, Pₛ)
solution = solve(problem, DPRKN6(), callback = cb, reltol=1e-6)

plot(solution, idxs=(2), linewidth=2,  framestyle = :box, yformatter = :scientific, xaxis = "t [s]", yaxis = "R(t) [m]")

David_F · October 3, 2022, 6:01am

Regarding my previous post, I now understand that there was no issue with updated values for Pᵥ and Pᵧ. After some light modifications, I think the code is OK. However, I encounter problems with different solvers. None of the stiff solvers work, all fail with the error “no method matching vec”, while ROCK4 and DPRKN6() work without problems. Tsit4 also produces the same problem.

Here’s the stack trace:

MethodError: no method matching vec(::Float64)
Closest candidates are:
  vec(::StrideArraysCore.PtrArray{S, D, T, 1, C, 0}) where {S, D, T, C} at C:\Users\DAVIDFUR\.julia\packages\StrideArraysCore\VQxXL\src\reshape.jl:1
  vec(::StrideArraysCore.PtrArray{S, D, T, N, C, 0}) where {S, D, T, N, C} at C:\Users\DAVIDFUR\.julia\packages\StrideArraysCore\VQxXL\src\reshape.jl:2
  vec(::StrideArraysCore.StrideArray) at C:\Users\DAVIDFUR\.julia\packages\StrideArraysCore\VQxXL\src\reshape.jl:11
  ...

Stacktrace:
  [1] copyto!(dest::Vector{Float64}, A::ArrayPartition{Float64, Tuple{Float64, Float64}})
    @ RecursiveArrayTools C:\Users\DAVIDFUR\.julia\packages\RecursiveArrayTools\IgoeN\src\array_partition.jl:171
  [2] copyto_axcheck!(dest::Vector{Float64}, src::ArrayPartition{Float64, Tuple{Float64, Float64}})
    @ Base .\abstractarray.jl:1127
  [3] Vector{Float64}(x::ArrayPartition{Float64, Tuple{Float64, Float64}})
    @ Base .\array.jl:626
  [4] (Vector)(x::ArrayPartition{Float64, Tuple{Float64, Float64}})
    @ Core .\boot.jl:481
  [5] zeromatrix(A::ArrayPartition{Float64, Tuple{Float64, Float64}})
    @ RecursiveArrayTools C:\Users\DAVIDFUR\.julia\packages\RecursiveArrayTools\IgoeN\src\array_partition.jl:371
  [6] build_J_W(alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, u::ArrayPartition{Float64, Tuple{Float64, Float64}}, uprev::ArrayPartition{Float64, Tuple{Float64, Float64}}, p::Function, t::Float64, dt::Float64, f::DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, #unused#::Type{Float64}, #unused#::Val{false})
    @ OrdinaryDiffEq C:\Users\DAVIDFUR\.julia\packages\OrdinaryDiffEq\vfMzV\src\derivative_utils.jl:889
  [7] alg_cache(alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, u::ArrayPartition{Float64, Tuple{Float64, Float64}}, rate_prototype::ArrayPartition{Float64, Tuple{Float64, Float64}}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, uprev::ArrayPartition{Float64, Tuple{Float64, Float64}}, uprev2::ArrayPartition{Float64, Tuple{Float64, Float64}}, f::DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, t::Float64, dt::Float64, reltol::Float64, p::var"#7#8", calck::Bool, #unused#::Val{false})
    @ OrdinaryDiffEq C:\Users\DAVIDFUR\.julia\packages\OrdinaryDiffEq\vfMzV\src\caches\rosenbrock_caches.jl:506
  [8] __init(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#7#8", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::ContinuousCallback{typeof(condition), typeof(affect!), typeof(affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Float64, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Rational{Int64}, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, 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, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq C:\Users\DAVIDFUR\.julia\packages\OrdinaryDiffEq\vfMzV\src\solve.jl:316
  [9] __solve(::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#7#8", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, ::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :reltol, :maxiters), Tuple{ContinuousCallback{typeof(condition), typeof(affect!), typeof(affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, Float64, Int64}}})
    @ OrdinaryDiffEq C:\Users\DAVIDFUR\.julia\packages\OrdinaryDiffEq\vfMzV\src\solve.jl:5
 [10] solve_call(_prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#7#8", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; merge_callbacks::Bool, kwargshandle::KeywordArgError, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :reltol, :maxiters), Tuple{ContinuousCallback{typeof(condition), typeof(affect!), typeof(affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, Float64, Int64}}})
    @ DiffEqBase C:\Users\DAVIDFUR\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:472
 [11] solve_up(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#7#8", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, sensealg::Nothing, u0::ArrayPartition{Float64, Tuple{Float64, Float64}}, p::Function, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :reltol, :maxiters), Tuple{ContinuousCallback{typeof(condition), typeof(affect!), typeof(affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, Float64, Int64}}})
    @ DiffEqBase C:\Users\DAVIDFUR\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:834
 [12] solve(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#7#8", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :reltol, :maxiters), Tuple{ContinuousCallback{typeof(condition), typeof(affect!), typeof(affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, Float64, Int64}}})
    @ DiffEqBase C:\Users\DAVIDFUR\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:801
 [13] top-level scope
    @ In[4]:82
 [14] eval
    @ .\boot.jl:368 [inlined]
 [15] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:1428

And here’s the updated code, which works with ROCKS4, but not for the other mentioned solvers:

using DifferentialEquations
using Plots                                          
                                                                                                                                                                            
#Parameters                                                     
Y = 2.0e7                       #Yield stress (N/m2)                       
µ = 4.54e9                       #shear modulus (N/m2)                     
ρ = 1840.0                      #density (Kg.m-3)                 
η = 1.0e2                       #dynamic viscosity ((N/m2).s)                 
                                                                                                            
#Porosity
φ = 0.01

#Initial conditions                                                              
a₀ = 10e-6                      # bubble inner radius (m)
b₀ = a₀/cbrt(φ)                 # bubble outer radius (m)
da₀ = 1e-7                      # initial radial velocity (ms-1)            
tspan = (0.0, 5.0e-6)           # time (s)         

#Distention ratios
α₀ = b₀^3/(b₀^3 - a₀^3)
α₁ = (2*µ*α₀ + Y)/(2*µ+Y)
α₂ = 2*µ*α₀/(2*µ+Y)
#print("α₀ = ",α₀,"  α₁=",α₁,"  α₂=",α₂)

Pₛ = t->1.0e16*t

# function definitions
function α(a, b)
   b^3/(b^3 - a^3)   
end

function b_radius(a)
   cbrt(b₀^3 - a₀^3 + a^3)
end

function c_radius(α)
   cbrt((2*µ*(a₀^3)*(α₀-α))/(Y*(α₀-1)))
end

function Pᵥ(a, c, da, new_α)
    if (new_α <= α₀) && (new_α >= α₁)
        0
    elseif (new_α <= α₁) && (new_α >= α₂)
        12*a*a*da*η*(-(1/3)*((1/c^3 - 1/a^3)))
    elseif (new_α <= α₂) && (new_α >= 1)
        12*a*a*da*η*(-(1/3)*((1/b_radius(a)^3) - 1/a^3))        
    else
        0
    end
end

function Pᵧ(a, new_α)
    if (new_α <= α₀) && (new_α >= α₁)
        4*µ*(α₀-new_α)/(3*new_α*(new_α-1))
    elseif (new_α <= α₁) && (new_α >= α₂)
        (2*Y/3)*((1-(2*µ)/(Y*new_α) + log((2*µ*(α₀-new_α))/(Y*(new_α-1)))))
    elseif (new_α <= α₂) && (new_α >= 1)
        (2*Y/3)*log(new_α/(new_α-1))
    else
        0
    end
end

function rpnnp(da, a, p, t)
    new_α = α(a, b_radius(a))
    #@show a, b_radius(a), new_α
    c = c_radius(α(a, b_radius(a)))
    @. (((-p(t)+Pᵧ(a, new_α)-Pᵥ(a, c, da, new_α))/ρ)-0.5*da*(4*da*(1-cbrt(φ))-da*(1-cbrt(φ)^4)))/(a*(1-cbrt(φ)))
end

function condition(y, t, integrator)
   #rise time of shock
   t >= 1.0e-16
end

function affect!(integrator)
   integrator.p[1] = 1.0e9
end
cb = ContinuousCallback(condition, affect!)

problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, Pₛ)
solution = solve(problem, Rodas4(), callback = cb, reltol=1e-15, maxiters=Int(1e6))

plot(solution, idxs=(2), linewidth=2,  framestyle = :box, yformatter = :scientific, xaxis = "t [s]", yaxis = "R(t) [m]", xlims = (0, 5), ylims = (0, 10))

I wonder if this could be a solver(s) issue?

David_F · October 18, 2022, 10:55am

I managed to solve the issue. Apparently, I had to define the scalars as vectors for the other solvers to work.

Now, I’m facing another issue - the callbacks to change the parameter seem to not work. I don’t receive any errors. The solver solves it. However, the callbacks aren’t changing the parameter.

I’d appreciate help with this.

The updated code is below:

using DifferentialEquations
using Plots
using MTH229

#Parameters                                                     
Y = 2.0e8                        #Yield stress (N/m2)                       
µ = 4.54e9                       #shear modulus (N/m2)                     
ρ = 1840.0                       #density (Kg.m-3)                 
η = 1.0e2                        #dynamic viscosity ((N/m2).s)                 
                                                                                                            
#Porosity
φ = 0.01

#Initial conditions                                                              
a₀ = [10e-6]                     #bubble inner radius (m)
b₀ = a₀./cbrt(φ)                 #bubble outer radius (m)
da₀ = [1e-4]                     #initial radial velocity (m*s⁻¹)            
tspan = (0.0, 10.0e-6)           #time (s)         

#Distention ratios
α₀ = b₀.^3 ./(b₀.^3 - a₀.^3)
α₁ = (2*µ.*α₀ .+ Y)./(2*µ.+Y)
α₂ = 2*µ.*α₀/(2*µ.+Y)

# function definitions
function α(a, b)
   b.^3 ./(b.^3 .- a.^3)
end

function b_radius(a)
   cbrt.(b₀.^3 .- a₀.^3 .+ a.^3)
end

function c_radius(α)
   cbrt.((2*µ.*(a₀.^3)./Y).*(α₀.-α)./(α₀.-1))
end

function Pᵥ(a, c, da, α)
    f(r) = 1/r^4
    @syms r
    if (α[1] <= α₀[1]) && (α[1] >= α₁[1])
        0
    elseif (α[1] <= α₁[1]) && (α[1] >= α₂[1])
        #12*a[1]^2*da[1]*η*(-(1/3)*((1/c[1]^3 - 1/a[1]^3)))       
        12*a[1]^2*da[1]*η*riemann(f, a[1], c[1], 100, method="simpsons")
    elseif (α[1] <= α₂[1]) && (α[1] >= 1)
        #12*a[1]^2*da[1]*η*(-(1/3)*((1/b_radius(a)[1]^3) - 1/a[1]^3))    
        12*a[1]^2*da[1]*η*riemann(f, a[1], b_radius(a)[1], 100, method="simpsons")
    else
        0
    end
end

function Pᵧ(α)
    if (α[1] <= α₀[1]) && (α[1] >= α₁[1])
        4*µ*(α₀[1]-α[1])/(3*α[1]*(α[1]-1.0))
    elseif (α[1] <= α₁[1]) && (α[1] >= α₂[1])
        2*Y/3*(1-(2*µ*(α₀[1]-α[1]))/(Y*α[1]) + log( abs((2*µ*(α₀[1]-α[1]))/(Y*(α[1]-1))) ))
    elseif (α[1] <= α₂[1]) && (α[1] >= 1)
        2*Y/3*log(α[1]/(α[1]-1))
    else
        0
    end
end

function rpnnp(dda, da, a, p, t)
    @. dda = ((-p(t) + Pᵧ(α(a, b_radius(a))) - Pᵥ(a, c_radius(α(a, b_radius(a))), da, α(a, b_radius(a))))/(ρ) - 0.5*da*(4*da*(1-cbrt(φ))-da*(1-cbrt(φ)^4)))/(a*(1-cbrt(φ)))
end

P = t->1.0e16*t   #external wave Pa*s⁻¹*s

function condition(y, t, integrator)
   #rise time of shock 0.1µs
   t >= 1.0e-7
end

function affect!(integrator)
   #constant external drive value Pa
   integrator.p[1] = 1.0e9
end
cb = ContinuousCallback(condition, affect!)

problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, P)
solution = solve(problem, Rodas4(), callback = cb, reltol = 1.0e-10)


plot(solution, xlims = (0, 4.0e-6), ylims = (0, 1.0e-5), idxs=(2), linewidth=2,  framestyle = :box, yformatter = :scientific, xaxis = "t [s]", yaxis = "R(t) [m]")

David_F · October 18, 2022, 12:39pm

I bypassed the problem by using a piecewise function instead of the callbacks.
Here’s the new piece of code:

P = t-> t<1.0e-7 ? 1.0e16*t : 1.0e9   #external wave profile Pa*s⁻¹*s

and the callback functions were removed.