PDE discretization fails when a parameter is used

[RickDW]

Hi all, I've working on a 2D hydrogen flame model that I've described in #59. I've tried to change the diffusivity parameter kappa by using @parameters x y t kappa and defining the PDESystem using

@named pdesys = PDESystem(
    eqs, # partial differential equations
    bcs, # initial/boundary conditions
    domains, # domain of the independent variables (i.e. spatial/time)
    [x,y,t], # independent variables
    [YF(x,y,t), YO(x,y,t), YP(x,y,t), T(x,y,t)], # dependent variables
    [κ], # parameters
    ) 

The error that I get is this:

Discretization failed at structural_simplify, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

ERROR: ArgumentError: Dict(kv): kv needs to be an iterator of tuples or pairs
Stacktrace:
 [1] Dict(kv::Vector{Any})
   @ Base ./dict.jl:132
 [2] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:124
 [3] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [4] top-level scope
   @ ~/Museum/uni/msc5/UQDDM/project/scripts/hydrogenFlame2.jl:142

caused by: BoundsError: attempt to access Num at index [2]
Stacktrace:
 [1] indexed_iterate(I::Num, i::Int64, state::Nothing)
   @ Base ./tuple.jl:95
 [2] grow_to!(dest::Dict{Num, Float64}, itr::Vector{Any}, st::Int64)
   @ Base ./dict.jl:153
 [3] grow_to!(dest::Dict{Any, Any}, itr::Vector{Any})
   @ Base ./dict.jl:145
 [4] dict_with_eltype
   @ ./abstractdict.jl:539 [inlined]
 [5] Dict(kv::Vector{Any})
   @ Base ./dict.jl:129
 [6] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:124
 [7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [8] top-level scope
   @ ~/Museum/uni/msc5/UQDDM/project/scripts/hydrogenFlame2.jl:142

[RickDW]

And here's the full script for completeness:

using ModelingToolkit
using MethodOfLines
using DomainSets
using DifferentialEquations
using GLMakie

## Define variables and derivative operators

# the two spatial dimensions and time
@parameters x y t κ
# the mass fractions MF (respectively H₂, O₂, and H₂O), temperature,
# and the source terms of the mass fractions and temperature
@variables YF(..) YO(..) YP(..) T(..) 
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

# laplace operator
∇²(u) = Dxx(u) + Dyy(u)

## Define domains

xmin = ymin = 0.0 # cm
xmax = 1.8 # cm
ymax = 0.9 # cm
tmin = 0.0
tmax = 0.06 # s

domains = [
    x ∈ Interval(xmin, xmax),
    y ∈ Interval(ymin, ymax),
    t ∈ Interval(tmin, tmax)
]

## Define parameters

# diffusivity coefficient (for temperature and mass fractions)
# κ =  2.0 #* cm^2/s
# constant (divergence-free) velocity field
U = [50.0, 0.0] #.* cm/s
# y coordinates of the inlet area on the left side of the domain
inlety = (0.3, 0.6) # mm

## Define the model

# diffusion and convection terms for the mass fractions and temperature
eqs = [
    Dt( YF(x,y,t) ) ~ κ * ∇²( YF(x,y,t) ) - U[1] * Dx( YF(x,y,t) ) - U[2] * Dy( YF(x,y,t) )
    Dt( YO(x,y,t) ) ~ κ * ∇²( YO(x,y,t) ) - U[1] * Dx( YO(x,y,t) ) - U[2] * Dy( YO(x,y,t) )
    Dt( YP(x,y,t) ) ~ κ * ∇²( YP(x,y,t) ) - U[1] * Dx( YP(x,y,t) ) - U[2] * Dy( YP(x,y,t) )
    Dt( T(x,y,t) )  ~ κ * ∇²( T(x,y,t) )  - U[1] * Dx( T(x,y,t) )  - U[2] * Dy( T(x,y,t) )
]

function atInlet(x,y,t)
    return (inlety[1] < y) * (y < inlety[2])
end

bcs = [
    #
    # initial conditions
    #

    T(x, y, 0) ~ 300.0, # K; around 26 C
    # domain is empty at the start
    YF(x,y,0) ~ 0.0,
    YO(x,y,0) ~ 0.0, 
    YP(x,y,0) ~ 0.0,

    #
    # boundary conditions
    #

    # left side
    T(xmin,y,t) ~ atInlet(xmin,y,t) * 950 + (1-atInlet(xmin,y,t)) *  300, # K
    YF(xmin,y,t) ~ atInlet(xmin,y,t) * 0.0282, # mass fraction
    YO(xmin,y,t) ~ atInlet(xmin,y,t) * 0.2259,
    YP(xmin,y,t) ~ 0.0,

    # right side
    Dx( T(xmax,y,t) ) ~ 0.0, # K/s
    Dx( YF(xmax,y,t) ) ~ 0.0, # mass fraction/s
    Dx( YO(xmax,y,t) ) ~ 0.0,
    Dx( YP(xmax,y,t) ) ~ 0.0,

    # top
    Dy( T(x,ymax,t) ) ~ 0.0, # K/s
    Dy( YF(x,ymax,t) ) ~ 0.0, # mass fraction/s
    Dy( YO(x,ymax,t) ) ~ 0.0,
    Dy( YP(x,ymax,t) ) ~ 0.0,

    # bottom
    Dy( T(x,ymin,t) ) ~ 0.0, # K/s
    Dy( YF(x,ymin,t) ) ~ 0.0, # mass fraction/s
    Dy( YO(x,ymin,t) ) ~ 0.0,
    Dy( YP(x,ymin,t) ) ~ 0.0,
]

@named pdesys = PDESystem(
    eqs, # partial differential equations
    bcs, # initial/boundary conditions
    domains, # domain of the independent variables (i.e. spatial/time)
    [x,y,t], # independent variables
    [YF(x,y,t), YO(x,y,t), YP(x,y,t), T(x,y,t)], # dependent variables
    [κ], # parameters
) 

## Discretize the system

N = 3 
dx = (xmax-xmin)/N
dy = (ymax-ymin)/N
order = 2

discretization = MOLFiniteDifference(
    [x=>dx, y=>dy], t, approx_order=order
)

# this creates an ODEProblem or a NonlinearProblem, depending on the system
problem = discretize(pdesys, discretization)

[xtalax]

Parameters need default values, like

@named pdesys = PDESystem(
    eqs, # partial differential equations
    bcs, # initial/boundary conditions
    domains, # domain of the independent variables (i.e. spatial/time)
    [x,y,t], # independent variables
    [YF(x,y,t), YO(x,y,t), YP(x,y,t), T(x,y,t)], # dependent variables
    [κ => 2.0], # parameters
) 

[RickDW]

Okay, thank you. Parameters need to be fully specified before discretization then it seems, is that right?

[xtalax]

At present yes, in fact one of my next tasks was to work out how to do parameter optimization with MethodOfLines and support that/write a tutorial, so you are a bit ahead of me in that respect. It may already be possible, you can get the ODESystem with symbolic_discretize(pdesys, discretization), in which case the parameters don't need default values.

[xtalax]

you could try defining a function of kappa including the discretize, solve and value extraction, then use GalacticOptim.jl

[RickDW]

Alright, I'll give that a go tomorrow. I don't need to optimize anything luckily, just need a few thousand evaluations at different parameter values. By the way, is there a way to access the system of ODEs of a discretized PDESystem? I realized today that I'll need access to that to create a surrogate model. If I understand you correctly, you get that for free if you use symbolic_discretize, but I haven't found a way to do that for the normal `discretize'. (by this I mean that I want to access the ODE system in the form of some kind of equations. My surrogate modeling approach requires splitting the discretized equations into a linear and a non-linear part)

[xtalax]

You can get the equations with

sys = symbolic_discretize(pdesys, disc)
simplesys = structural_simplify(sys)
eqs = equations(simplesys)
@show eqs

if you look at https://github.com/SciML/MethodOfLines.jl/blob/1e099d65ba473f86300723725f9b61be42b77206/src/discretization/MOL_discretization.jl#L151 you can see that's pretty much exactly what we're doing when generating the problem. Will you generate the surrogate automatically from the ModelingToolkit equations? If so please share your code, we want that as a method.

[xtalax]

I was incorrect before, the problem does not need to be re-discretized every time. You can do

prob_with_param(kappa) = remake(prob, p = [κ => kappa])
newprob = prob_with_param(rand())
solve(newprob, ...)

[RickDW]

Ah I remember seeing the remake function before in the docs, that's a nice catch.

If I can get the automatic surrogate generation up and running I'll happily share it with you. I imagine it could be useful to many different people. I haven't fully wrapped my head around the theory yet, and I'm on a tight deadline, but I'll let you know what happens. In case you're curious, the paper I'm reproducing uses Proper Orthogonal Decomposition and the Discrete Empirical Interpolation Method. It's called "Projection-based model reduction for reacting flows" by Buffoni and Willcox.

Just now I tried to run symbolic_discretize without passing default parameter values to the system constructor, and it's thrown this error again:

Discretization failed at structural_simplify, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

ERROR: ArgumentError: Dict(kv): kv needs to be an iterator of tuples or pairs
Stacktrace:
 [1] Dict(kv::Vector{Any})
   @ Base ./dict.jl:132
 [2] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:119
 [3] top-level scope
   @ ~/Museum/uni/msc5/UQDDM/project/scripts/hydrogenFlameSteady.jl:131

caused by: BoundsError: attempt to access Num at index [2]
Stacktrace:
 [1] indexed_iterate(I::Num, i::Int64, state::Nothing)
   @ Base ./tuple.jl:95
 [2] grow_to!(dest::Dict{Any, Any}, itr::Vector{Any})
   @ Base ./dict.jl:142
 [3] dict_with_eltype
   @ ./abstractdict.jl:539 [inlined]
 [4] Dict(kv::Vector{Any})
   @ Base ./dict.jl:129
 [5] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:119
 [6] top-level scope
   @ ~/Museum/uni/msc5/UQDDM/project/scripts/hydrogenFlameSteady.jl:131

Maybe I misunderstood you, but I don't think this is supposed to happen with symbolic_discretize?

[xtalax]

I just noticed that the error message is wrong, it was failing at the system construction step. You do have to supply default parameter values but they are not substituted in to the equations at the symbolic_discretize stage so you can still do this manually with the eqs.

Could you share the paper?

[RickDW]

Ah okay, that makes a lot of sense. I'm looking right at the equations resulting from symbolic_discretize, which are still working with the unassigned parameter.

The paper is freely available in google scholar, that might be easier to find than sharing a pdf myself.

[RickDW]

Small remark about discretize's implementation: structural_simplify is applied only when symbolic_discretize outputs an ODESystem and not when it's a NonlinearProblem. For my model I get 100 equations out of symbolic_discretize but structural_simplify manages to reduce it to 36. I think it would be more efficient to be working with the simplified system at solve time (?)

[ChrisRackauckas]

yes we should apply it in both cases.

[xtalax]

@RickDW Were you able to get this working? How did it go?

[RickDW]

@xtalax Hey, things are going relatively smoothly. I've got the code that reduces the order of a nonlinear set of equations set up, the next step is to hook up this code so I can automatically generate a new nonlinear problem from the equations. I think modelingtoolkit supplies an interface for that, so that might be relatively easy to code up.

What is probably going to take the most time is translating between solutions of the reduced order model and the original high order model. I haven't looked into that yet though, so I'll see what the best approach is. Once this is implemented the surrogate generation code is pretty much fully functional, although it might have some rough edges. It feels like this might be similar to the solution extraction code that we talked about before. You mentioned that there's some work being done on that, so maybe I can help out.

My plans have changed a bit though, meaning that I'll probably take a break from implementing these final few bits for the next two weeks or so. It's become quite a big project so I'll pace myself once I start again. The implementation is my first priority, with proper documentation right after that. Once things are stable we can start thinking about integrating it into a sciml package. Speaking of which, would you want this included in modeling toolkit? It could also be put in surrogates.jl, but since my implementation only supports symbolic systems (it could be extended to non-symbolic systems, but I don't think I'll do that myself) modeling toolkit might be a better fit. What do you think about this?

(Right now my code is not public yet, but I'll post it in this repository soon)

[xtalax]

@ChrisRackauckas Where would this ^ go?

[ChrisRackauckas]

Try MTK. I've been wanting forms of model order reduction there (there's a whole GSoC about it). We'll see if it fits.

[RickDW]

Alright then, I'll let you know once I've got something to integrate.

[xtalax]

@RickDW How did you get remake to work? Having trouble getting this to work for an example

[RickDW]

I think the problem for me was the way the parameters were passed to remake.

I was incorrect before, the problem does not need to be re-discretized every time. You can do

prob_with_param(kappa) = remake(prob, p = [κ => kappa])
newprob = prob_with_param(rand())
solve(newprob, ...)

Instead of using a dictionary like you did here, I'm using a vector. It's kind of annoying since you need to know the order of the parameters but at least it works.