Problems with MethodOfLines discretizing PDESystem involving equations of symbolic arrays

General Usage
, , ,
1

Hi, I’m using ModelingToolkit to do some work involving Riemannian geometry, and in particular would like to work with tensor fields in arbitrary dimension. The following is as close to a MWE as I can get my code, I think:

using Symbolics, MethodOfLines, LinearAlgebra, DomainSets, ModelingToolkit
dim = 3
@parameters p[1:dim]
@variables u(..) (g(..))[1:dim, 1:dim]
Dp = [Differential(pi) for pi in p]
# some boundary conditions for u on the boundary of a cube
bcs = [
    u(0, [p[i] for i in 2:dim]...) ~ 0,
    u(1, [p[i] for i in 2:dim]...) ~ 1,
    [
        u([p[j] for j in 1:(i - 1)]..., 0, [p[j] for j in (i + 1):dim]...) ~ p[1]
        for i in 2:dim
    ]...,
    [
        u([p[j] for j in 1:(i - 1)]..., 1, [p[j] for j in (i + 1):dim]...) ~ p[1]
        for i in 2:dim
    ]...,
]
domains = [(pi ∈ Interval(0, 1)) for i in 1:dim]

eqs = [
    g(p...) .~ I(dim) * (1 + 0.5 * sin(p[1])),
    # just a Laplacian equation, the exact form is not important, 
    # but it involves the metric tensor g,
    # which may be expensive to evaluate
    sum([
        (Dp[i](inv(g(p...))[i, j] * sqrt(abs(det(g(p...)))) * Dp[j](u(p...)))) for
        i in 1:dim, j in 1:dim
    ]) ~ 0,
]
@named PDEsys = PDESystem(eqs, bcs, domains, p, [u(p...)])

dpi = 0.01
discretizer = MOLFiniteDifference([pi => dpi for pi in p])
prob = discretize(PDEsys, discretizer)

Creating the PDESystem works fine, but discretizing seems to get stuck on the first of the two eqs, namely the one involving the auxiliary variable g. It throws the error:

ERROR: LoadError: type Arr has no field lhs
Stacktrace:
  [1] getproperty(x::Symbolics.Arr{Any, 2}, f::Symbol)
    @ Base ./Base.jl:37
  [2] split_complex_eq(eq::Symbolics.Arr{Any, 2}, redvmaps::Vector{Pair{…}}, imdvmaps::Vector{Pair{…}})
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/make_pdesys_compatible.jl:58
  [3] (::PDEBase.var"#132#140")(eq::Symbolics.Arr{Any, 2})
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/make_pdesys_compatible.jl:109
  [4] _mapreduce(f::PDEBase.var"#132#140", op::typeof(vcat), ::IndexLinear, A::Vector{Any})
    @ Base ./reduce.jl:440
  [5] _mapreduce_dim(f::Function, op::Function, ::Base._InitialValue, A::Vector{Any}, ::Colon)
    @ Base ./reducedim.jl:365
  [6] mapreduce(f::Function, op::Function, A::Vector{Any})
    @ Base ./reducedim.jl:357
  [7] handle_complex(pdesys::PDESystem)
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/make_pdesys_compatible.jl:108
  [8] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/symbolic_discretize.jl:9
  [9] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, kwargs::@Kwargs{})
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/discretization_state.jl:58
 [10] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/discretization_state.jl:55

Stepping through with the Debugger seems to show that checking whether the equation is complex with !isequal(real(eq),eq) in hascomplex (further up in handle_complex returns true for the broadcasted equation. If I instead do not broadcast the equation involving g (i.e. use ~ instead of .~), I get an error

ERROR: LoadError: MethodError: no method matching real(::SymbolicUtils.BasicSymbolic{Matrix{Real}})
Stacktrace:
  [1] hascomplex(term::SymbolicUtils.BasicSymbolic{Matrix{Real}})
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/symbolic_utils.jl:292
  [2] hascomplex(eq::Equation)
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/symbolic_utils.jl:291
  [3] #127
    @ ~/.julia/packages/PDEBase/imf8O/src/make_pdesys_compatible.jl:86 [inlined]
  [4] _any(f::PDEBase.var"#127#135", itr::Vector{Equation}, ::Colon)
    @ Base ./reduce.jl:1220
  [5] any(f::Function, a::Vector{Equation}; dims::Function)
    @ Base ./reducedim.jl:1020
  [6] handle_complex(pdesys::PDESystem)
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/make_pdesys_compatible.jl:86
  [7] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/symbolic_discretize.jl:9
  [8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, kwargs::@Kwargs{})
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/discretization_state.jl:58
  [9] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
    @ PDEBase ~/.julia/packages/PDEBase/imf8O/src/discretization_state.jl:55

in that same hascomplex function. Am I missing something here? I’m still unsure if there is some way to tell PDESystem that I want to use g as an auxiliary variable, maybe that would fix something? Just inserting the expression for g in the differential equation for u directly (and removing the equation for g) seems to work, but I am summing over possibly quite a few dimensions and each evaluation of g may be expensive, so I’d like to avoid that.

2

I noticed that what I’m doing here is probably not how one is supposed to do this. I guess, I should not declare g as its own unknown variable but just do something like @register_array_simbolic g(p...) = ..., and then let the discretization automatically remove any duplicate computations of g. Can anyone confirm if this is the right approach? Should I maybe even do the same to avoid lots of computations of inv(g) occurring in the sum defining the second equation?

Apart from that I think I’ll just close this for now.

3

It shouldn’t be required, but right now not all optimization passes exist so doing some manually can help.

1 Like