One-dimensional NLS using MethodOfLines ‘type Array has no field lhs`

[ruvilecamwasam]

I’m trying to solve the one-dimensional nonlinear Schrodinger equation using MethodOfLines.jl. I looked at the code on the docs: http://methodoflines.sciml.ai/dev/tutorials/brusselator/ 3, and adapted it as:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters x z
@variables A(..)
Dx   = Differential(x)
Dz   = Differential(z)
Dzz  = Differential(z)^2

xmin = 0.
xmax = 1e-1
zmax = 10.
zmin = -zmax

c0 = 1.
A0(x,z) = c0*sech(c0*z/sqrt(2))*exp(im*c0^2*x/2)

domains = [x ∈ Interval(xmin,xmax), z ∈ Interval(zmin,zmax)]

eq = [im*Dx(A(x,z))+Dzz(A(x,z)) ~ -abs2(A(x,z))*A(x,z)]
bcs = [A(xmin,z) ~ A0(xmin,z), 
        A(x,zmin) ~ 0, 
        A(x,zmax) ~ 0]

@named pdesys = PDESystem(eq,bcs,domains,[x,z],[A(x,z)])

N       = 100
dz      = 1/N
order   = 2

discretization = MOLFiniteDifference([z=>dz], x)

@time prob = discretize(pdesys,discretization)

However this won’t run, giving error type Array has no field lhs. I think the only main difference between my code and the tutorial is that my PDE is for a single function.

I'm running Julia v1.7.2, with package versions:

• DiffEqOperators v4.41.0
• MethodOfLines v0.2.0
• ModelingToolkit v6.7.1 ModelingToolkit v8.5.5
• DomainSets v0.5.9

The full error is:

ERROR: type Array has no field lhs
Stacktrace:
  [1] getproperty
    @ .\Base.jl:42 [inlined]
  [2] (::MethodOfLines.var"#65#67"{Vector{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}}})(x::Vector{Equation})
    @ MethodOfLines C:\Users\ruvil\.julia\packages\MethodOfLines\9MLPx\src\MOL_utils.jl:64
  [3] mapreduce_first(f::MethodOfLines.var"#65#67"{Vector{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}}}, op::Function, x::Vector{Equation})
    @ Base .\reduce.jl:394
  [4] _mapreduce(f::MethodOfLines.var"#65#67"{Vector{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}}}, op::typeof(union), #unused#::IndexLinear, A::Vector{Vector{Equation}})
    @ Base .\reduce.jl:405
  [5] _mapreduce_dim(f::Function, op::Function, #unused#::Base._InitialValue, A::Vector{Vector{Equation}}, #unused#::Colon)
    @ Base .\reducedim.jl:330
  [6] #mapreduce#725
    @ .\reducedim.jl:322 [inlined]
  [7] mapreduce(f::Function, op::Function, A::Vector{Vector{Equation}})
    @ Base .\reducedim.jl:322
  [8] get_all_depvars(pdesys::PDESystem, depvar_ops::Vector{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}})
    @ MethodOfLines C:\Users\ruvil\.julia\packages\MethodOfLines\9MLPx\src\MOL_utils.jl:64
  [9] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
    @ MethodOfLines C:\Users\ruvil\.julia\packages\MethodOfLines\9MLPx\src\discretization\MOL_discretization.jl:18
 [10] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
    @ MethodOfLines C:\Users\ruvil\.julia\packages\MethodOfLines\9MLPx\src\discretization\MOL_discretization.jl:146
 [11] top-level scope
    @ .\timing.jl:220

[ChrisRackauckas]

How are you using those versions? MOL v0.2 requires ModelingToolkit v8. It won't let you install that set.

[ruvilecamwasam]

Sorry, I mistranscribed. I am running ModelingToolkit v8.5.5 (I double checked and the other version numbers I wrote are correct).

[xtalax]

Complex equations are not currently supported but I will be working on this, can you reformulate as a system of equations in terms of the real and imaginary parts of
A?

[ruvilecamwasam]

Thanks, I gave it a try but it doesn't seem to be working.

The PDE is the nonlinear Schrödinger equation:

im Dx(A)+Dzz(A)=-|A|^2A,

where Dx=d/dx, and Dzz=d^2/dz^2. Defining A(x,z)=a(x,z)+im b(x,z) for real functions a(x,z) and b(x,z), this becomes

-Dx(b)+Dzz(a)=-(a^2+b^2)a
Dx(a)+Dzz(b) = -(a^2+b^2)b

Following the tutorial for the Brusselator, I coded this up as below (discretising over z):

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters x z
@variables a(..) b(..)
Dx   = Differential(x)
Dz   = Differential(z)
Dzz  = Differential(z)^2

xmin = 0.
xmax = 1e-1
zmax = 10.
zmin = -zmax

c0 = 1.
A0(x,z) = c0*sech(c0*z/sqrt(2))*exp(im*c0^2*x/2) #Initial condition for complex PDE
A0r(x,z) = c0*sech(c0*z/sqrt(2))*cos(c0^2*x/2) #Real part of initial condition
A0i(x,z) = c0*sech(c0*z/sqrt(2))*sin(c0^2*x/2) #Imaginary part of initial condition

domains = [x ∈ Interval(xmin,xmax), z ∈ Interval(zmin,zmax)]

eq = [-Dx(b(x,z))+Dzz(a(x,z)) ~ -(a(x,z)^2+b(x,z)^2)*a(x,z),
       Dx(a(x,z))+Dzz(b(x,z)) ~ -(a(x,z)^2+b(x,z)^2)*b(x,z)]

bcs = [a(xmin,z) ~ A0r(xmin,z), 
        b(xmin,z) ~ A0i(xmin,z),
        a(x,zmin) ~ 0, b(x,zmin) ~ 0,
        a(x,zmax) ~ 0, b(x,zmax) ~ 0]

@named pdesys = PDESystem(eq,bcs,domains,[x,z],[a(x,z),b(x,z)])

N       = 10 #Also tried running for N=100, got same result
dz      = 1/N
order   = 2

discretization = MOLFiniteDifference([z=>dz], x, approx_order=order)

@time prob = discretize(pdesys,discretization)

This gave me a very long error, they key parts of are below. I can't seem to find the problem in my code though.

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

ERROR: ExtraVariablesSystemException: The system is unbalanced. There are 402 highest order derivative variables and 398 equations.
More variables than equations, here are the potential extra variable(s):

Stacktrace:
 [1] error_reporting(state::TearingState{ODESystem}, bad_idxs::Vector{Int64}, n_highest_vars::Int64, iseqs::Bool)
   @ ModelingToolkit.StructuralTransformations C:\Users\ruvil\.julia\packages\ModelingToolkit\a7NhS\src\structural_transformation\utils.jl:35
 [2] check_consistency(state::TearingState{ODESystem})
   @ ModelingToolkit.StructuralTransformations C:\Users\ruvil\.julia\packages\ModelingToolkit\a7NhS\src\structural_transformation\utils.jl:66
 [3] structural_simplify(sys::ODESystem; simplify::Bool)
   @ ModelingToolkit C:\Users\ruvil\.julia\packages\ModelingToolkit\a7NhS\src\systems\abstractsystem.jl:923
 [4] structural_simplify(sys::ODESystem)
   @ ModelingToolkit C:\Users\ruvil\.julia\packages\ModelingToolkit\a7NhS\src\systems\abstractsystem.jl:920
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines C:\Users\ruvil\.julia\packages\MethodOfLines\9MLPx\src\discretization\MOL_discretization.jl:151
 [6] top-level scope
   @ .\timing.jl:220

[xtalax]

This seems to be above board, this is a bug. I'll take a closer look at this soon.

[ruvilecamwasam]

Thanks! Let me know if you want me to try anything else.