Modifying the input parameters to function in a function while solving differential equations

New to Julia
1

I am new to Julia. Some workshops from this year’s JuliaCon really impressed me and I started exploring… This might be quite a basic question.

I am trying to implement a callback while solving an ODE problem. I want the callback to be based on a value of the passed parameter which I want to modify in the function itself. This is my function.

function df(du, u, (bool_burst, bp, p), t)
    .....
    
    press = #some calculated value
 
    bool_burst = press<bp
    
    
    .....
    
end


function condition(u, t, integrator)
    @show integrator.p[1]
    integrator.p[1]==1
    
end 

function affect!(integrator)
    .....
end

In my case, the value of bool_burst is not changing when I use @show to print the value.
Thanks in advance.

2

It is a little hard to tell without seeing more code, but with what you posted, I do not expect the change in bool_burst to be visible outside of du. I think you need to add press to your state vector u and then use a continuous callback like

function continuous_condition(u, t, integrator)
   press = u[some index]
   return press - integrator.p[1]
end

Then affect! will be called when continuous_condition crosses zero.

1 Like
3

Thank You for your response @contradict.
Actually press is standard atmospheric pressure calculated using a separate function, at given altitude which is tracked in du.

function df(du, u, (bp,p), t)
    x = u[1]
    y = u[2]
    z = u[3]
    vx = u[4]
    vy = u[5]
    vz = u[6]
    Tib = u[7]
    Tf = u[8]
    mg = u[9]
    
    a = atm(floor(y))
    
    ## here press is calculated
    ρ , press, T, ν = ISAdata(y)

    ####
    
    g = 9.81
    
    
    
    rn = 3* mg *p["mol_m_a"]
    rd = 4*π*ρ*p["mol_m_g"]
    r = (rn/rd)^(1.0/3)
    area = π*r^2
    
    Re = ρ*sqrt(vx^2 + vy^2 + vz^2)*2*r/ν
    cib = press/(8314.0 * Tib)
    ρg = cib*p["mol_m_g"]
    
    Sca = ν/(ρ*2*r)
    Scg = Sca
    
    
    acc_y = get_forces_x( p["mp"] + mg, mg, ρ, g, p["cd_up"], vy, area, p["mol_m_g"], p["mol_m_a"])/(mg+ p["mp"])
    dm = get_dm(T, press, 2*r, Re, Scg, Sca, cib, p)
    dTf, dTib = DT(Tf, T, Tib, Re, 2*r, ρg, g, 0.0, press, y, t, mg, p)

    du[1] = 0
    du[2] = vy
    du[3] = 0
    du[4] = 0
    du[5] = acc_y
    du[6] = 0
    du[7] = dTib
    du[8] = dTf
    du[9] = -dm
    
end


function check(u, t, integrator)
   ##Here we have to calculate pressure (press) again to compare, but it is extra calculation for no reason
    press = atm(u[2]).pressure #same as ISAdata used above but this gives only pressure.
    p =  press - integrator.p[1]
    p
    
end

I wanted to not calculate press again inside check .

Your answer will definitely work. If there is any other way where I don’t need to track press in du. That would be great.

Thank you.

1 Like
4

I see, I think recompute is the way to go here. In the presence of adaptive time stepping, attempting to keep track of the computed value independent of the state vector would be quite difficult. I would not be surprised to learn that the solver might call check() at a time for which an interpolated value of u was used, so no corresponding call to df() existed.

1 Like
5

Thank you @contradict. It was very helpful. I was not aware of the fact that check() can be called even without df(), which seems quite possible.
Thank you.