# Solve ode problem with parameters described as function of the variables

Hello all. I have a sinusoidal excitational model of the hydraulic pumps.
I have an expression which calculates the density, viscosity and bulk modulus as a fluid(P,T) ,P meaning pressure is a variable.

``````function fluid(Pre::Float64,Temp::Float64)
χ₀=0.001;                                 # Air content percentage in the oil at STP
Pₛₐₜ=4e6;                                  # Saturation pressure of the oil at stp
Pᵥₐₚ=1e4;                                 # Vapor pressure of the oil at stp (Shell's data)
αₜ=7.2e-4;                                # Volumetric thermal expansion coeff.  (/K)
ρₗ₀=872.0;                                   # oil density at stp
ρₐ₀=1.225;                                  # air density at stp
νₚₐ=60.3095;           # Kinematic viscosity of liquid oil at operating temperature.(cSt) (Solved by params.m file)

bulk_liq=(10^((0.3766*(log10(νₚₐ)^0.3307))-0.2766)
+(5.851-0.01382*Temp)*(Pre/1e9))*1e9;
μₐ=(1.458e-6*(Temp+273.15)^1.5)/(Temp+110.4+273.15);        #Sutherland Expression for air
μₗ₀=(νₚₐ*1e-6)*ρₗ₀*exp(-αₜ*(Temp-T_atm::Float64));              #Dynamic viscosity of liquid oil
if Pre>=Pₛₐₜ
rho_f=(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)*exp((Pre-p_atm::Float64 )/bulk_liq-αₜ*(Temp-T_atm::Float64));
Beta_f=bulk_liq;
mu_f=(νₚₐ*1e-6)*(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)*exp((Pre-p_atm::Float64 )/bulk_liq-αₜ*(Temp-T_atm::Float64));
elseif Pre<Pₛₐₜ && Pre>Pᵥₐₚ
phi=(Pₛₐₜ-Pre)/(Pₛₐₜ-Pᵥₐₚ);
rho_f=(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)/
((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/273.15)*(p_atm::Float64 /Pre)^(1/1.42))+exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)));

#Beta_f=(1+(χ₀)/(1-χ₀)*((Temp+273.15)/(273.15))*((p_atm::Float64 /Pre)^(1/1.42))*phi)/
#(bulk_liq^-1 + (χ₀)/(1-χ₀)*((Temp+273.15)/(273.15))*(p_atm::Float64 /Pre)^(1/1.42)*(1/(Pₛₐₜ-p_atm::Float64 ))*((Pₛₐₜ-Pre)/1.42/Pre+1) );

Beta_f=((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/(273.15))*((p_atm::Float64 /Pre)^(1/1.42)))+exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)))/
(((χ₀)/(1-χ₀)*((Temp+273.15)/(T_atm::Float64 +273.15))*(p_atm::Float64 /Pre)^(1/1.42))*((1/(Pₛₐₜ-Pᵥₐₚ))+(phi/(1.42*Pre)))+
((exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)))/bulk_liq));

mu_f=((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/273.15)*((p_atm::Float64 /Pre)^(1/1.42))*μₐ)+
(exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64))*μₗ₀))/((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/273.15)*((p_atm::Float64 /Pre)^(1/1.42)))+
exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)));
elseif Pre<=Pᵥₐₚ
rho_f=(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)/
((χ₀/(1-χ₀)*((Temp+273.15)/273.15)*((p_atm::Float64 /Pᵥₐₚ)^(1/1.42)))+exp((p_atm::Float64 -Pᵥₐₚ)/bulk_liq+αₜ*(Temp-T_atm::Float64)));
Beta_f=1.42*Pre;
mu_f=μₐ;
end
return rho_f,mu_f,Beta_f
end
``````

Creating my ode function, I did this:

``````function damper!(dPdt,P,pars,t)
dₚ,dᵣ,Aₚ,Aᵣ,Aₜₚ,hₚ,hᵣ,Δₚ,Δᵣ,lₚₛ,δₚₛ,lᵣₛ,δᵣₛ,
Hₜ,
dₒco,dₒre,dₒba,lₒco,lₒre,lₒba,Cfmax_co,Cfmax_re,Cfmax_ba,
rᵥₐco,lᵥₐco,rₑₓₜco,Aₑₓₜco,kᵥₐco,x₀co,xₛₜᵣco,lₛₚₒco,δco,
rᵥₐre,lᵥₐre,rₑₓₜre,Aₑₓₜre,kᵥₐre,x₀re,xₛₜᵣre,lₛₚₒre,δre,
rᵥₐba,lᵥₐba,rₑₓₜba,Aₑₓₜba,kᵥₐba,x₀ba,xₛₜᵣba,lₛₚₒba,δba,
dₚₐ,lₚₐ,Cfmaxₕ,d_f₁,d_f₂,A_bot,Aₜₒₚ,Cfₗₐₙ,kᵢₙ,x₀,
V₀co,V₀re,Vg₀,Pg₀,Vres₀,
amp,freq,=pars;

#--Pressure in the 3 Chambers--#
Pᵣₑ,Pcₒ,Pᵣₑₛ=P;                     # Pressure in the rebound chamber,compression chamber and reservoir

#--Imposed Sinusodial Excitation-#
x=amp*sin(2*pi*freq*t);             # Displacement from the middle of the piston
u=amp*2*pi*freq*cos(2*pi*freq*t);   # Velocity of the piston rod

ρ₁,μ₁,β₁=fluid(Pᵣₑ,T);
ρ₂,μ₂,β₂=fluid(Pcₒ,T);
ρ₃,μ₃,β₃=fluid(Pᵣₑₛ,T);

Q_12=ori(Pᵣₑ,Pcₒ,dₒco,lₒco,Cfmax_co,ρ₁,μ₁)+
ori(Pᵣₑ,Pcₒ,dₒre,lₒre,Cfmax_re,ρ₁,μ₁)+
sum(spring.(Pᵣₑ,Pcₒ,ρ₁,rᵥₐre,lᵥₐre,rₑₓₜre,Aₑₓₜre,kᵥₐre,x₀re,xₛₜᵣre,lₛₚₒre,δre,3))+
pi*dₚ*δₚₛ^3/12/lₚₛ/μ₁*abs(Pᵣₑ-Pcₒ)*(Pᵣₑ>=Pcₒ);              #Leakage of the piston seals

Q_21=ori(Pcₒ,Pᵣₑ,dₒco,lₒco,Cfmax_co,ρ₂,μ₂)+
ori(Pcₒ,Pᵣₑ,dₒre,lₒre,Cfmax_re,ρ₂,μ₂)+
sum(spring.(Pcₒ,Pᵣₑ,ρ₂,rᵥₐco,lᵥₐco,rₑₓₜco,Aₑₓₜco,kᵥₐco,x₀co,xₛₜᵣco,lₛₚₒco,δco,3))+
pi*dₚ*δₚₛ^3/12/lₚₛ/μ₂*abs(Pcₒ-Pᵣₑ)*(Pcₒ>=Pᵣₑ);             #Leakage of the piston seals

Q_23=ori(Pcₒ,Pᵣₑₛ,dₒba,lₒba,Cfmax_ba,ρ₂,μ₂)+
sum(spring.(Pcₒ,Pᵣₑₛ,ρ₂,rᵥₐba,lᵥₐba,rₑₓₜba,Aₑₓₜba,kᵥₐba,x₀ba,xₛₜᵣba,lₛₚₒba,δba,3));

Q_32=ori(Pᵣₑₛ,Pcₒ,dₒba,lₒba,Cfmax_ba,ρ₃,μ₃)+
check(Pᵣₑₛ,Pcₒ,Pg₀,dₚₐ,lₚₐ,Cfmaxₕ,d_f₁,d_f₂,A_bot,Aₜₒₚ,Cfₗₐₙ,kᵢₙ,x₀);

dPdt[1]=β₁/(V₀re-Aₜₚ*x)*(Q_21*(ρ₂/ρ₁)-Q_12
+Aₜₚ*u);

dPdt[2]=β₂/(V₀co+Aₚ*x)*(Q_12*(ρ₁/ρ₂)-Q_21+Q_32*(ρ₃/ρ₂)-Q_23
-Aₚ*u);

dPdt[3]=(Q_23*(ρ₂/ρ₃)-Q_32)/
(((Vres₀+Vg₀-((Pg₀/Pᵣₑₛ)^(1/1.42))*Vg₀)/β₃)+(Pg₀^(1/1.42)*Vg₀/1.42/Pᵣₑₛ^(1+1/1.42)));

nothing

#Force=P_ex*A_tp-P_co*A_p+pi*vis_dyna((P_co+P_ex)/2)*(d_p*h_p/Delta_p+0.5*d_r*h_gu/Delta_r)*abs(u);
end
``````

While solving the ode with Rodas5(), I get an error with a huge page, few lines of it being:

``````MethodError: no method matching fluid(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}, ::Float64)
Closest candidates are:
fluid(::Float64, ::Float64) at In[89]:1

Stacktrace:
[1] damper!(dPdt::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}}, P::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}}, pars::Vector{Any}, t::Float64)
@ Main .\In[94]:20
[2] (::ODEFunction{true, SciMLBase.AutoSpecialize, typeof(damper!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing})(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}}, ::Vararg{Any})
``````

I guess I called the fluid function in a matlab way which would give real values of density, viscosity and bulk modulus. I have checked coupled of tutorials with this kind of a problem where a parameter is function of the variable but I didn’t get a grasp.

Your help will be highly appreciated.

Welcome!

In Julia you don’t need to specify argument types in the method signature.

So in this case I would suggest the following signature

``````function fluid(Pre, Temp)
``````

what’s happening here is the solver by default tries to use auto differentiation of your code which fails because of the Float64s (which in your case don’t do anything). for models that actually can’t be autodiffed, you can pass autodiff=false to the solver, but here is better to just fix your model.

I read couple of posts on how to optimise a function, and one said I should specify the input variables type.
I tried to remove all of the `::Float64`, but still the same problem.

It seems yes, but still I don’t understand how does that affect a scalar input when the variable had an explicit value before performing the differentiation.
I tried to put Rodas5(autodiff=false) and it seems to take ages.

The next problem you’re having is that `T_atm` and `p_atm` are global variables. Pass them in to `fluid` and everything will be a lot faster. I’m pretty sure you mis-interpreted the advice about specifying types.

Since the error is specifically linked to the method signature, I don’t believe that is correct. You might need to restart julia to override the original method definition (it will take precendence if you don’t).

After performing the calculations with p_atm and T_atm as local variables in the `fluid(P)`.

``````Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).
└ @ SciMLBase C:\Users\490132\.julia\packages\SciMLBase\MVhUm\src\integrator_interface.jl:502
``````

I tried CVODE_BDF, Rodas4 and Vern7 got the same answer.

Also, what matlab solver were you using? It’s possible there’s a better Julia solver for your problem than `Rodas5` (`Rodas5P` is almost always better and depending on tolerances and a bunch of other stuff you might want to be using `QNDF`/`FBDF` or `Rosenbrok23`).

Is there a reason you are using stiff solvers? I would start with `Tsit5()` if you don’t need a stiff solver

the fact that OP is hitting maxiters implies that the problem is likely stiff.

I used ode23s to start with. It didn’t give me great performance.

The problem is as stiff as a hard bread would be early morning. I can’t tell, I tried with `abstol=1e-8` and `reltol=1e-6`.

2 Likes

What did you use in Matlab?

Method to solve the ode? ode23s and RKF45. The time it takes to calculate especially this check valve equation;

``````function check(P_res,P_co,Pg_0,d_o,l_open,Cfmax_o,d_film_1,d_film_2,A_bot,A_top,C_lang,k_in,h_0)
rho,vis=fluid(P_res/2+P_co/2)[1:2];
if P_res*A_bot-P_co*A_top-Pg_0*(A_top-A_bot)-k_in*h_0>=0
Cdv_max=0.27;                                        #Annular Flow maximum value.
f(h::Float64) = P_res*A_bot - P_co*A_top - k_in*(h+h_0)+
2*(P_res-P_co)*(C_lang/(pi*d_o^2/4)-A_bot/2*(Cfmax_o*pi*d_o^2/4)^-2)*
((Cfmax_o*pi*d_o^2/4)^-2+(Cdv_max*tanh(4*h/vis/100*sqrt(2*(P_res-P_co)*rho))*pi*d_film_1*h)^-2)^-1;
h_sol=find_zero(f,0.0001,Order2());
Q_valve=sqrt(2*(P_res-P_co)/rho*
((Cfmax_o*pi*d_o^2/4)^-2+(Cdv_max*tanh(4*h_sol/vis/100*sqrt(2*(P_res-P_co)*rho))*pi*d_film_2*h_sol)^-2)^-1);
else
Q_valve=0;
end
end
``````

I am sick and tired of this modeling problem.

Rosenbrok23 is the Julia equivalent of ode23s, so it’s probably a decent place to start. How do the Julia and matlab solve times compare? It should be pretty easy to get a fairly substantial speedup while keeping the solvers the same from the language change, but that will require a bit of work to improve your Julia code.

Tried Rosenbrok23 too, same problem with maxiters. I think the problem is the functions that am calling in the ode by `ori(),spring(),check()` where it should have its way of declaring them.
I am no expert, I have searched everywhere on the internet and I just found one example of parameter as a function of time not a variable.

I recommend going through How to debug Julia simulation codes (ODEs, optimization, etc.!) - YouTube to see if any of the debugging techniques there help. There’s a pretty high likelyhood you’ve swapped a sign somewhere that makes the model unstable.

1 Like

Have watched this JuliaCon and used almost everything being said. I think the problem is the discontinuity by the exogenous functions being called in the damper function. And I can’t make them smooth as they should be.

Can you predict at which points the exogenous functions are not smooth?