Solve ode problem with parameters described as function of the variables

New to Julia
, ,
1

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.

2

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)
3

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.

4

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.

5

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.

6

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.

7

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).

8

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

I have received this error about Larger maxiters,

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.

9

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).

10

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

11

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

12

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

13

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
14

What did you use in Matlab?

15

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. :ok_man:

16

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.

17

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.

18

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
19

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.

20

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