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

Specific Domains Modelling & Simulations
, , ,
24

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