Dependent mass matrix with RadauIIA5

General Usage
1

When running this code Im obtaing not expect not smooth solutions.

  using Plots
  using DifferentialEquations

  Ta=20.0+273.15;                                                               
  ka=1.5207e-11*Ta^3 - 4.8574e-8*Ta^2 + 1.0184e-4*Ta - 3.9333e-4;                                     
 ro_a=360.77819/Ta; 
 eta_a=(-1.1555e-14*Ta^3 + 9.5728e-11*Ta^2 + 3.7604e-8*Ta - > 3.4484e-6)*ro_a;   
  Q=0.6/3600;                
  R0=0.00135;                
  L=0.140;                   
  T0=150.0+273.15;            
  ro=781.5;                        
  N0=3.0;                     
  eta_T0=2.2e4;              
  lambda_T0=2.0;             
  Ea=5.77e4;                
  alpha=0.105;               
  g=9.81;                     
  u0=Q/(ro*3.1415*R0^2);   
  G0=eta_T0/lambda_T0; 
  V1=0.42*L*ka/(ro*u0^(2/3)*R0^(5/3))*(ro_a/eta_a)^(1/3); 
  V3=G0/(ro*T0); 
  C1=ro*u0^2/G0;
  C3=(ro*g*L)/G0; 
  Eeq=(1.0+N0)/(2.0*N0)+((9.0+6.0*N0+9.0*N0^2.0)^0.5/6.0/N0);
  czz0=2.0;

function filament(du,u,p,t)

    Cp=1980+3.7*T0*u[5] 
    EE=(9.0*N0-2.0*u[1]-u[2])/(9.0*N0-6.0*u[1]-3.0*u[2])
    dEE=18.0*N0/(9.0*N0-6.0*u[1]-3.0*u[2])^2
    lambda=lambda_T0*exp(Ea/(8.314*T0)*((1.0/u[5])-1.0))
    Cf=0.205*((ro_a*u[4]*u0*u[3]*R0)/eta_a)^(-0.61)
    C2=2.0*Cf*ro_a*L*u0^2/(G0*R0)   

    println(u[3])

    du[1] = u[5]/(lambda*u0/L)*((1.0-alpha) + alpha*EE*u[1]/u[5])*(EE*u[1]/u[5]-1.0)
    du[2] = u[5]/(lambda*u0/L)*((1.0-alpha) + alpha*EE*u[2]/u[5])*(EE*u[2]/u[5]-1.0)
    du[3] = 0.0
    du[4] = -C2*u[4]^2/u[3] + C3
    du[5] = -V1/(Cp*u[4]^(2/3)*u[3]^(5/3))*(u[5]-Ta/T0)

    nothing
end    

function update_func(M,u,p,t)

    Cp=1980+3.7*T0*u[5] 
    EE=(9.0*N0-2.0*u[1]-u[2])/(9.0*N0-6.0*u[1]-3.0*u[2])
    dEE=18.0*N0/(9.0*N0-6.0*u[1]-3.0*u[2])^2
    lambda=lambda_T0*exp(Ea/(8.314*T0)*((1.0/u[5])-1.0)) 
    Cf=0.205*((ro_a*u[4]*u0*u[3]*R0)/eta_a)^(-0.61) 
    C2=2.0*Cf*ro_a*L*u0^2/(G0*R0)

    M = [-u[4]                      0                       0                   -u[1]                       0
            0                     -u[4]                     0               +2.0*u[2]                       0
            0                       0                     -u[4]               -u[3]/2.0                       0
        -2.0*dEE*(u[2]-u[1])+EE   -(dEE*(u[2]-u[1])+EE)   -2.0*EE*(u[2]-u[1])/u[3]   +C1*u[4]                       0

            0                       0                       0            -V3*EE*(u[2]-u[1])/(Cp*u[4])      +1.0]

end

    uâ‚€ = [1.0/Eeq, czz0, 1.0, 1.0, 1.0] 
    tspan = (0.0,1.0)
    
    Cp0=1980+3.7*T0*uâ‚€[5]      #969.9 
    EE0=(9.0*N0-2.0*uâ‚€[1]-uâ‚€[2])/(9.0*N0-6.0*uâ‚€[1]-3.0*uâ‚€[2]) 
    dEE0=18.0*N0/(9.0*N0-6.0*uâ‚€[1]-3.0*uâ‚€[2])^2 
    lambda0=lambda_T0*exp(Ea/(8.314*T0)*((1.0/uâ‚€[5])-1.0)) 
    Cf0=0.205*((ro_a*uâ‚€[4]*u0*uâ‚€[3]*R0)/eta_a)^(-0.61) 
    C20=2*Cf0*ro_a*L*u0^2/(G0*R0)      

    M0 = [-uâ‚€[4]                      0                         0                           -uâ‚€[1]                       0 
                0                     -uâ‚€[4]                     0                          +2.0*uâ‚€[2]                       0 
                0                       0                     -uâ‚€[4]                        -uâ‚€[3]/2                       0 
        -2.0*dEE0*(uâ‚€[2]-uâ‚€[1])+EE0   -(dEE0*(uâ‚€[2]-uâ‚€[1])+EE0)   -2.0*EE0*(uâ‚€[2]-uâ‚€[1])/uâ‚€[3]   +C1*uâ‚€[4]                       0 
                0                       0                       0                   -V3*EE0*(uâ‚€[2]-uâ‚€[1])/(Cp0*uâ‚€[4])      +1]

    M = DiffEqArrayOperator(M0,update_func=update_func)

    f = ODEFunction(filament,mass_matrix=M)
    prob_mm = ODEProblem(f,uâ‚€,tspan)

    sol = DifferentialEquations.solve(prob_mm, RadauIIA5())

    display(plot(sol, tspan=(0.0, 1.0), layout=(5,1)))

Besides, there is a message

Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(filament),DiffEqArrayOperator{Float64,Array{Float64,2},typeof(update_func)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64}}(1.0,0.0,0.0,0.0,1.0,0.0)Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(filament),DiffEqArrayOperator{Float64,Array{Float64,2},typeof(update_func)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64}}(1.0,0.0,0.0,1.0,0.0,0.0)

Is something wrong here?

2

Dependent mass matrices like that are pretty fundamentally incompatible with automatic differentiation (if you want to do it in a way that’s fast), so after experimenting a lot with this, we are going to be dropping the experimental support for dependent mass matrices. That said, any dependent mass matrix problem can be rewritten as a constant mass matrix problem, which is a numerically much better approach (and something we hope to automate through ModelingToolkit fairly soon).

1 Like
3

Thanks Chris.
If I understood correctly, converting the system to the form du=M-1(t,u)f(t,u) will give better reuslts.
In addition I ve rewritten the system as f(t,u,du)=0 and runned Sundials IDA; the result is good for some conditions. However, depending on the value of parameters the solution can diverge very quickly.

4

No, don’t invert it (that won’t be possible with DAEs anyways). Just write the version that’s the equivalent equation with constant mass matrices through substitution.

5

Got it! I going to try that way.
Thanks