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.