# Problems with MethodOfLines discretizing PDESystem involving equations of symbolic arrays
**URL:** https://discourse.julialang.org/t/problems-with-methodoflines-discretizing-pdesystem-involving-equations-of-symbolic-arrays/119916
**Category:** General Usage
**Tags:** question, package, modelingtoolkit, methodoflines
**Created:** [September 26, 2024, 12:29pm UTC](https://discourse.julialang.org/t/problems-with-methodoflines-discretizing-pdesystem-involving-equations-of-symbolic-arrays/119916 "2024-09-26T12:29:08Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@r0uv3n](https://discourse.julialang.org/u/r0uv3n)
#### Post date: [September 26, 2024, 12:29pm UTC](https://discourse.julialang.org/t/problems-with-methodoflines-discretizing-pdesystem-involving-equations-of-symbolic-arrays/119916/1 "2024-09-26T12:29:09Z")
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:
```julia
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:
```julia
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
```julia
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.
---
### Author:  [@r0uv3n](https://discourse.julialang.org/u/r0uv3n)
#### Post date: [September 28, 2024, 9:07am UTC](https://discourse.julialang.org/t/problems-with-methodoflines-discretizing-pdesystem-involving-equations-of-symbolic-arrays/119916/2 "2024-09-28T09:07:49Z")
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 28, 2024, 10:14am UTC](https://discourse.julialang.org/t/problems-with-methodoflines-discretizing-pdesystem-involving-equations-of-symbolic-arrays/119916/3 "2024-09-28T10:14:21Z")
> [@r0uv3n](#):
>
> and then let the discretization automatically remove any duplicate computations of `g`. Can anyone confirm if this is the right approach?
It shouldn’t be required, but right now not all optimization passes exist so doing some manually can help.