Interpolated parameters into the PDESystems

Specific Domains Modelling & Simulations
,
1

Hey I am trying to solve a system of equations where my input parameters are altitude dependent values. These values should be supplied into system of equations after discretizing the equation. I have created an Interpolated function which has to be taken into the equations and has to solve. I made an example code here. Your help is appreciated.

using ModelingToolkit, DomainSets
using MethodOfLines, OrdinaryDiffEq
using Interpolations
using Symbolics: Num
using RuntimeGeneratedFunctions

Independent variables

@parameters t x
@variables u(..)

Correct Symbolic Definition: Callable symbolic parameter

@parameters D_param(..)

Differentials

Dt = Differential(t)
Dx = Differential(x)

=== INTERPOLATION SETUP ===

Grid

x_vals = 0.0:0.1:Float64(pi)
const XMAX = Float64(pi)

Parameter data (Diffusion coefficient)

D_vals = ones(length(x_vals)) * 0.1

Create Interpolation object

D_fun = LinearInterpolation(x_vals, D_vals, extrapolation_bc=Line())

Correct parameter dictionary input:

param_vals = Dict(D_param(x) => x → D_fun(x))

Goal: While discretization D_param(x) value should be taken from interpolated function while solving

PDE equation: Only includes the D(x) term

eqs = [
Dt(u(t,x)) ~ D_param(x) * Dx(Dx(u(t,x)))
]

Domains

tdom = Interval(0.0, 1.0)
xdom = Interval(0.0, XMAX)
domains = [t ∈ tdom, x ∈ xdom]

Boundary Conditions (BCs)

bcs = [
u(0.0, x) ~ sin(x), # IC at t=0
u(t, 0.0) ~ 1e-3, # BC at x=0
u(t, XMAX) ~ 0.0 # BC at x=π
]

Build PDESystem

@named pdesys = PDESystem(eqs, bcs, domains, [t,x], [u(t,x)], [D_param(x)];
defaults = param_vals );

discretization = MOLFiniteDifference([x => 10], t;approx_order = 2, use_ODAE = false)
sym_prob = symbolic_discretize(pdesys, discretization); ## I could see the discretized equations over 10 altitudes.
prob = discretize(pdesys, discretization);

prob has threw the following error:
The system of equations is:
Equation[Differential(t)((u(t))[2]) - (8.207015875029361(u(t))[1] - 16.414031750058722(u(t))[2] …

MethodError: Cannot convert an object of type var"#110#111" to an object of type Number
The function convert exists, but no method is defined for this combination of argument types.

Closest candidates are:
convert(::Type{Number}, !Matched::Unitful.Quantity)
@ Unitful ~/.julia/packages/Unitful/sxiHA/src/conversion.jl:166
convert(::Type{T}, !Matched::Union{Static.StaticBool{N}, Static.StaticFloat64{N}, Static.StaticInt{N}} where N) where T<:Number
@ Static ~/.julia/packages/Static/SeEGr/src/Static.jl:427
convert(::Type{T}, !Matched::Number) where T<:Number
@ Base number.jl:7
…

After this I have to do the following one
sol = solve(prob, Tsit5())

2

Interpolations.jl isn’t setup with the Symbolics.jl tooling yet. If you want to integrate an interpolation, you’ll need to use DataInterpolations.jl. That only does 1-dimensional interpolations though. If you need higher dimension interpolations then you’ll want GitHub - SciML/DataInterpolationsND.jl: Interpolation of arbitrarily high dimensional array data, though this integration is not complete yet:

When that PR is merged this story for ND will be complete.

3

Thanks for the quick help @ChrisRackauckas. I presented simple case of a large system. With Datainterpolation and I could solve the given example. But in the actual problem, I am solving equations using a reaction network coupled with method of lines. This requires Datainterpolation of input parameters multiple times depending on the chosen spatial grid that make my calculation very slow. Wondering if there is a better approach to avoid this Datainterpolation. Please let me know if any tutorial example. Thanks.

4

Are you asking for the ND interpolation?

5

thanks for the help. I understand that parameters have to be interpolated and the interpolated functions have to plugged into the equations system to solve the system of equations. I have constructed PDESystem of 15 equations with 84 parameters which are interpolated functions using DataInterpolations. But at the step prob = discretize(pdesys, discretization); it is taking unusual time for N = 45 spatial points. Wondering if there is a faster approach for this? Please suggest me.

6

At this time no, we’re working on that exact problem right now though

7

Ok. thanks…!!