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

**URL:** <https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777>\
**Category:** Modelling & Simulations\
**Tags:** performance, pde, modelingtoolkit, differentialequation\
**Created:** [February 19, 2022, 10:05pm UTC](https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777 "2022-02-19T22:05:44Z")\
**Posts on this page:** 7\
**Page:** 2

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 1, 2022, 12:08am UTC](https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777/21 "2022-03-01T00:08:59Z")

</div>

> [@RobertS](#):
>
> Unfortunately, I get an `UndefVarError: x201 not defined.` Why is that? What can I do to get it to work?

You didn’t supply initial conditions.

---

<div class="post-metadata">

**Author:** ![RobertS](https://avatars.discourse-cdn.com/v4/letter/r/ee7513/32.png) [@RobertS](https://discourse.julialang.org/u/RobertS)\
**Post date:** [March 1, 2022, 1:30pm UTC](https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777/22 "2022-03-01T13:30:40Z")

</div>

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/](https://mtk.sciml.ai/stable/mtkitize_tutorials/modelingtoolkitize/)

Nevertheless, I tried the following:

```julia
# 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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 1, 2022, 1:38pm UTC](https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777/23 "2022-03-01T13:38:14Z")

</div>

Can you share what the error is?

---

<div class="post-metadata">

**Author:** ![RobertS](https://avatars.discourse-cdn.com/v4/letter/r/ee7513/32.png) [@RobertS](https://discourse.julialang.org/u/RobertS)\
**Post date:** [March 1, 2022, 8:24pm UTC](https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777/24 "2022-03-01T20:24:15Z")

</div>

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

Complete code to reproduce:

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 6, 2022, 3:13am UTC](https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777/25 "2022-03-06T03:13:14Z")

</div>

Thanks, solution and tests here: [https://github.com/SciML/ModelingToolkit.jl/pull/1479](https://github.com/SciML/ModelingToolkit.jl/pull/1479)

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [March 30, 2022, 7:46pm UTC](https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777/26 "2022-03-30T19:46:48Z")

</div>

If simplicity of syntax and ease of use are what you are after, check out [MethodOfLines.jl](https://github.com/SciML/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.

---

<div class="post-metadata">

**Author:** ![RobertS](https://avatars.discourse-cdn.com/v4/letter/r/ee7513/32.png) [@RobertS](https://discourse.julialang.org/u/RobertS)\
**Post date:** [March 31, 2022, 7:30pm UTC](https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777/27 "2022-03-31T19:30:44Z")

</div>

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.

[Previous page](https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777.md?page=1)