Help me outperform Matlab in numerical solution of semidiscretized PDE system as much as possible

Don’t I unpack these variables from the parameter argument passed to the DAE function and hence, they are local? Or am I missing something?

Oh you’ve right, this should be fine

sparse(M). Do it if it’s big and sparse.

What’s more important is Jacobian sparsity. Did you read this tutorial?

https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/

1 Like

Thank you, I must have missed this part of the tutorial earlier. It seems promising.

I tried:

using Symbolics
du0 = copy(u0)
jac_sparsity = Symbolics.jacobian_sparsity((du,u)->DiscretizedSystemMassMatrixForm!(du,u,p,0.0),du0,u0)

This gives me the following error:

MethodError: no method matching Float64(::Num)

I suspect that this is due to how I defined parts of my parameter vector as dualcache. If I change my DAE function and parameter vector to not use dual cache variables at all, getting the sparsity works.

What can I do to get above lines of code working with my parameter vector definition and DAE function?

Just have an allocating version or something to make it easy.

Thanks, I followed your advice and indeed was able to derive and use the jacobian sparsity pattern for a decent speed up.

I was able to do so by having a second version of the DAE function which does not use dual cache, but otherwise is a 100% copy of the original DAE function. This is of course quite ugly, as I now have pretty much the same function twice. Is there a better way?

Additionally, please allow me to repeat some questions from my original post that have not been adressed in the thread yet. Would love to get some feedback on those :slight_smile:

It’s ugly, and we really just need to improve our sparsity interpreter to avoid it. @shashi can you think of any tricks that would work?

This would be the same issue as the sparsity one.

In general, the FAQ and PSA are helpful.

https://diffeq.sciml.ai/stable/basics/faq/

Thanks a lot. That is a lot of material to digest.

Thanks for the hint. I tried to do exactly that: I used a version of the DAE function which does not use dual caches and tried to modellingtoolkitize the ODEProblem defined by that DAE function:

using DifferentialEquations
using LinearAlgebra

# Generate the constants
N = 100
dx = 0.5 / N
q = dualcache(zeros(N,1))
rR = dualcache(zeros(N,1))
ϱu = 0
dp = 100e2
Tin = 200

# generate index lists to find the solution variables easily in solution vector
indices_Tgas = 1:N
indices_Tsol = indices_Tgas .+ N
indices_mflux = indices_Tsol[end] + 1
indices_RConcentrationSolid = indices_mflux[end] .+ (1:N)

indices = (indices_Tgas, indices_Tsol, indices_mflux, indices_RConcentrationSolid)

p = (N, # discretization intervals
     dx, # control volume length
     dp, # boundary condition: pressure differential across domain
     Tin, # boundary condition: Gas inlet temperature
     q, # caching variable for q
     rR, # caching variable for rR
     ϱu, # caching variable for ϱu (mass flux)
     indices_Tgas, # indices for gas temperature in solution vector
     indices_Tsol, # indices for solids temperature in solution vector
     indices_mflux,
     indices_RConcentrationSolid) # index for mass flux in solution vector


function DiscretizedSystemMassMatrixForm_noDualcache!(du,u,p,t)
  N,dx, dp, Tin, q, rR, ϱu, indices_Tgas, indices_Tsol, indices_mflux, indices_RConcentrationSolid = p
  Tgas = @view u[indices_Tgas]
  Tsol = @view u[indices_Tsol]
  dTgas = @view du[indices_Tgas]
  dTsol = @view du[indices_Tsol]
  mflux = @view u[indices_mflux]
  dmflux = @view du[indices_mflux]
  RConcsol  = @view u[indices_RConcentrationSolid]
  dRConcsol  = @view du[indices_RConcentrationSolid]
  
  q = 1e4 * (Tgas - Tsol)
  ϱu = mflux[1]

  rR = 5.5e6 .* exp.(.-5680 ./(273 .+ Tsol)) .* (1 .- exp.(.-RConcsol./50))

  # gas energy balance
  @inbounds begin
    dTgas[1] =  ϱu * 1000 * (Tin - Tgas[1]) - q[1] * dx
  end
  @inbounds for k in 2:N
    dTgas[k] = ϱu * 1000 * (Tgas[k-1] - Tgas[k]) - q[k] * dx
  end

  
  # solids energy balance
  @inbounds for k in 1:N
    dTsol[k] = q[k] / (1000 * 2000)
  end

  # gas momentum balance
  @inbounds begin
    dmflux[1] = - mflux[1] + 5e-3 * sqrt(dp)
  end

  # reactant concentration equations
  @inbounds for k in 1:N
    dRConcsol[k] = -rR[k]
  end
  
  
end

# build initial conditions
u0 = zeros(N*3 + 1)
u0[indices_Tgas] .= Tin
u0[indices_Tsol] .= 20
u0[indices_mflux] = 0.5
u0[indices_RConcentrationSolid] .= 3000

# build mass matrix
Ii = [indices_Tsol; indices_RConcentrationSolid]; 
Jj = [indices_Tsol; indices_RConcentrationSolid]; 
V = ones(length(Ii));
Msparse = sparse(Ii,Jj,V, 3*N+1, 3*N+1)


using ModelingToolkit
f_DAE = ODEFunction(DiscretizedSystemMassMatrixForm_noDualcache!, mass_matrix=Msparse)
probDAE = ODEProblem(f_DAE,u0,(0.0,1500),p)
de = modelingtoolkitize(probDAE)

Unfortunately, modelingtoolkitize still throws the same error as before:

ArgumentError: invalid index: α₈ of type Num

What am I doing wrong?

It looks like you can’t pass the index definitions through the parameters since ModelingToolkit wants to make all the parameters symbolic. I think this code is missing a couple of using statements too. This seems to work:

using DifferentialEquations
using DifferentialEquations.DiffEqBase: dualcache
using LinearAlgebra
using ModelingToolkit
using SparseArrays

# Generate the constants
const N = 100
dx = 0.5 / N
q = dualcache(zeros(N,1))
rR = dualcache(zeros(N,1))
ϱu = 0
dp = 100e2
Tin = 200

# generate index lists to find the solution variables easily in solution vector
const indices_Tgas = 1:N
const indices_Tsol = indices_Tgas .+ N
const indices_mflux = indices_Tsol[end] + 1
const indices_RConcentrationSolid = indices_mflux[end] .+ (1:N)

# indices = (indices_Tgas, indices_Tsol, indices_mflux, indices_RConcentrationSolid)

p = (# N, # discretization intervals
     dx, # control volume length
     dp, # boundary condition: pressure differential across domain
     Tin, # boundary condition: Gas inlet temperature
     q, # caching variable for q
     rR, # caching variable for rR
     ϱu, # caching variable for ϱu (mass flux)
     # indices_Tgas, # indices for gas temperature in solution vector
     # indices_Tsol, # indices for solids temperature in solution vector
     # indices_mflux,
     # indices_RConcentrationSolid, # index for mass flux in solution vector
)

function DiscretizedSystemMassMatrixForm_noDualcache!(du,u,p,t)
  # N,dx, dp, Tin, q, rR, ϱu, indices_Tgas, indices_Tsol, indices_mflux, indices_RConcentrationSolid = p
  dx, dp, Tin, q, rR, ϱu = p
  Tgas = @view u[indices_Tgas]
  Tsol = @view u[indices_Tsol]
  dTgas = @view du[indices_Tgas]
  dTsol = @view du[indices_Tsol]
  mflux = @view u[indices_mflux]
  dmflux = @view du[indices_mflux]
  RConcsol  = @view u[indices_RConcentrationSolid]
  dRConcsol  = @view du[indices_RConcentrationSolid]

  q = 1e4 * (Tgas - Tsol)
  ϱu = mflux[1]

  rR = 5.5e6 .* exp.(.-5680 ./(273 .+ Tsol)) .* (1 .- exp.(.-RConcsol./50))

  # gas energy balance
  @inbounds begin
    dTgas[1] =  ϱu * 1000 * (Tin - Tgas[1]) - q[1] * dx
  end
  @inbounds for k in 2:N
    dTgas[k] = ϱu * 1000 * (Tgas[k-1] - Tgas[k]) - q[k] * dx
  end


  # solids energy balance
  @inbounds for k in 1:N
    dTsol[k] = q[k] / (1000 * 2000)
  end

  # gas momentum balance
  @inbounds begin
    dmflux[1] = - mflux[1] + 5e-3 * sqrt(dp)
  end

  # reactant concentration equations
  @inbounds for k in 1:N
    dRConcsol[k] = -rR[k]
  end


end

# build initial conditions
u0 = zeros(N*3 + 1)
u0[indices_Tgas] .= Tin
u0[indices_Tsol] .= 20
u0[indices_mflux] = 0.5
u0[indices_RConcentrationSolid] .= 3000

# build mass matrix
Ii = [indices_Tsol; indices_RConcentrationSolid];
Jj = [indices_Tsol; indices_RConcentrationSolid];
V = ones(length(Ii));
Msparse = sparse(Ii,Jj,V, 3*N+1, 3*N+1)


f_DAE = ODEFunction(DiscretizedSystemMassMatrixForm_noDualcache!, mass_matrix=Msparse)
probDAE = ODEProblem(f_DAE,u0,(0.0,1500),p)
de = modelingtoolkitize(probDAE)
2 Likes

Thank you for the assistance. I got the MTKized system to work


# build initial conditions
u0 = zeros(N*3 + 1)
u0[indices_Tgas] .= Tin
u0[indices_Tsol] .= 20
u0[indices_mflux] = 0.5
u0[indices_RConcentrationSolid] .= 3000

# build mass matrix
Ii = [indices_Tsol; indices_RConcentrationSolid];
Jj = [indices_Tsol; indices_RConcentrationSolid];
V = ones(length(Ii));
Msparse = sparse(Ii,Jj,V, 3*N+1, 3*N+1)


f_DAE = ODEFunction(DiscretizedSystemMassMatrixForm_noDualcache!, mass_matrix=Msparse)
probDAE = ODEProblem(f_DAE,u0,(0.0,1500),p_noDualCache)
de = modelingtoolkitize(probDAE)

# MTKized system back to numerical problem:
probDAE_MTK1 = ODEProblem(de,[],(0.0,1500),p_noDualCache,jac=true,sparse=true)

# solve the numerical problem - works
sol_MTK1 = solve(probDAE_MTK1,QBDF())

# pretty much as fast as solving non-MTKized system with Jacobian sparsity, as in above posts.
@benchmark solve(probDAE_MTK1,QBDF()) 

To my surpise, the MTKized system does not solve faster than the original system with Jacobian sparsity information, although the MTKized system should use the analytic Jacobian, or am I mistaken?

Anyway. In a next step, I tried to structual_simplify the system and solve it:

# Now try to structural_simplify the system before solving it
de2 = structural_simplify(de)

probDAE_MTK2 = ODEProblem(de2,[],(0.0,1500),p_noDualCache,jac=true,sparse=true)

sol_MTK2 = solve(probDAE_MTK2,QBDF())

Unfortunately, I get an UndefVarError: x201 not defined. Why is that? What can I do to get it to work?

Thank you!

You didn’t supply initial conditions.

Are you sure this is the problem? As it seems, I don’t need to give the initial conditions explicity when defining the ODEProblem probDAE_MTK1 either. And this is in line with the following tutorial:
https://mtk.sciml.ai/stable/mtkitize_tutorials/modelingtoolkitize/

Nevertheless, I tried the following:

# Now try to structural_simplify the system before solving it
de2 = structural_simplify(de)

u0reduced = [u0[1:200]; u0[202:end]]
probDAE_MTK2 = ODEProblem(de2,u0reduced,(0.0,1500),p_noDualCache,jac=true,sparse=true)

sol_MTK2 = solve(probDAE_MTK2,QBDF())

This still gives the same error, which I find relieving, because defining u0reduced forced me to make assumptions about state / equation ordering in system de2 after structural_simplify, which can very well be wrong.

So, what is my problem?

Can you share what the error is?

I get the following error: UndefVarError: x_{201} not defined

Complete code to reproduce:

using DifferentialEquations
using DifferentialEquations.DiffEqBase: dualcache
using LinearAlgebra
using ModelingToolkit
using SparseArrays
using BenchmarkTools

# Generate the constants
const N = 100
dx = 0.5 / N
q = dualcache(zeros(N,1))
rR = dualcache(zeros(N,1))
ϱu = 0
dp = 100e2
Tin = 200

# generate index lists to find the solution variables easily in solution vector
const indices_Tgas = 1:N
const indices_Tsol = indices_Tgas .+ N
const indices_mflux = indices_Tsol[end] + 1
const indices_RConcentrationSolid = indices_mflux[end] .+ (1:N)

# indices = (indices_Tgas, indices_Tsol, indices_mflux, indices_RConcentrationSolid)

p_noDualCache = (dx, # control volume length
    dp, # boundary condition: pressure differential across domain
    Tin, # boundary condition: Gas inlet temperature
    q, # caching variable for q
    rR, # caching variable for rR
    ϱu)

function DiscretizedSystemMassMatrixForm_noDualcache!(du,u,p,t)
  # N,dx, dp, Tin, q, rR, ϱu, indices_Tgas, indices_Tsol, indices_mflux, indices_RConcentrationSolid = p
  dx, dp, Tin, q, rR, ϱu = p
  Tgas = @view u[indices_Tgas]
  Tsol = @view u[indices_Tsol]
  dTgas = @view du[indices_Tgas]
  dTsol = @view du[indices_Tsol]
  mflux = @view u[indices_mflux]
  dmflux = @view du[indices_mflux]
  RConcsol  = @view u[indices_RConcentrationSolid]
  dRConcsol  = @view du[indices_RConcentrationSolid]

  q = 1e4 * (Tgas - Tsol)
  ϱu = mflux[1]

  rR = 5.5e6 .* exp.(.-5680 ./(273 .+ Tsol)) .* (1 .- exp.(.-RConcsol./50))

  # gas energy balance
  @inbounds begin
    dTgas[1] =  ϱu * 1000 * (Tin - Tgas[1]) - q[1] * dx
  end
  @inbounds for k in 2:N
    dTgas[k] = ϱu * 1000 * (Tgas[k-1] - Tgas[k]) - q[k] * dx
  end


  # solids energy balance
  @inbounds for k in 1:N
    dTsol[k] = q[k] / (1000 * 2000)
  end

  # gas momentum balance
  @inbounds begin
    dmflux[1] = - mflux[1] + 5e-3 * sqrt(dp)
  end

  # reactant concentration equations
  @inbounds for k in 1:N
    dRConcsol[k] = -rR[k]
  end


end



# build initial conditions
u0 = zeros(N*3 + 1)
u0[indices_Tgas] .= Tin
u0[indices_Tsol] .= 20
u0[indices_mflux] = 0.5
u0[indices_RConcentrationSolid] .= 3000

# build mass matrix
Ii = [indices_Tsol; indices_RConcentrationSolid];
Jj = [indices_Tsol; indices_RConcentrationSolid];
V = ones(length(Ii));
Msparse = sparse(Ii,Jj,V, 3*N+1, 3*N+1)


f_DAE = ODEFunction(DiscretizedSystemMassMatrixForm_noDualcache!, mass_matrix=Msparse)
probDAE = ODEProblem(f_DAE,u0,(0.0,1500),p_noDualCache)
de = modelingtoolkitize(probDAE)

# MTKized system back to numerical problem:
probDAE_MTK1 = ODEProblem(de,[],(0.0,1500),p_noDualCache,jac=true,sparse=true)

# solve the numerical problem - works
sol_MTK1 = solve(probDAE_MTK1,QBDF())

# pretty much as fast as solving non-MTKized system with Jacobian sparsity, as in above posts.
@benchmark solve(probDAE_MTK1,QBDF()) 




# Now try to structural_simplify the system before solving it
de2 = structural_simplify(de)

# system de2 only has 200 states after structural simplification
probDAE_MTK2 = ODEProblem(de2,[],(0.0,1500),p_noDualCache,jac=true,sparse=true)

# Error thrown: UndefVarError: x_{201} not defined 
sol_MTK2 = solve(probDAE_MTK2,QBDF())

Thanks, solution and tests here: Invalidate system cache by ChrisRackauckas · Pull Request #1479 · SciML/ModelingToolkit.jl · GitHub

1 Like

If simplicity of syntax and ease of use are what you are after, check out MethodOfLines.jl, with which you can specify, discretize and solve PDE systems in close to purely mathematical syntax.

Compile times are quite high right now, but the solves themselves are comparable to manual semidiscretization, so still useful for parameter sweeps, ensembles and optimization problems.

Changes are coming which will improve compile times by a lot.

Thank you, I am more or less aware of the beauty the usage of MethodOfLines and ModelingToolkit offers, and I will definitely make use of them, when (if) I make the decision to use Julia in professional and productive settings. I am currently evaluating this.
For that reason, I wanted to get a feeling how much the performance of above codes can be pushed with some effort in coding.