Unclear error messages on BCs?

[henry2004y]

Hi,

I have been testing a PDE system with this package. This PDE system contains mixed differentials $$\frac{\partial^2 A_y}{\partial x\partial z}$$, which is currently not supported, as mentioned in #137. One suggested workaround is to add auxiliary variables. Here is what I have for the test problem:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

const Lx = 25.6
const Lz = 12.8
const nx = 8
const nz = 8
# background field
const B₀ = 1.0
# mass density
const ρ₀ = 1.0
# mass density at infinity
const ρ∞ = 0.2*ρ₀
# width of current sheet
const λ = 0.5
# perturbation amplitude of the magnetic flux
const ψ₀ = 0.1
# initial plasma β
const β = 1.0
# pressure normalization parameter
const p₀ = 0.5*β*B₀^2

# physical parameters
const γ = 5/3
const η = 0.0
const ν = 0.0


@parameters x z t
@variables ρ(..) p(..) ux(..) uz(..) Ay(..) Bx(..) Bz(..)
Dt = Differential(t)
Dx = Differential(x)
Dz = Differential(z)
Dxx = Differential(x)^2
Dzz = Differential(z)^2

∇²(u) = Dxx(u) + Dzz(u)

x_min = -Lx/2
z_min = -Lz/2
t_min = 0.0
x_max = Lx/2
z_max = Lz/2
t_max = 10.0

dx = Lx / nx
dz = Lx / nz

ψ(x,z,t) = ψ₀*cos(2π*x/Lx)*cos(π*z/Lz)

ρ0(x,z,t) = ρ₀*sech(z/λ)^2 + ρ∞

p0(x,z,t) = begin
   b = B₀*tanh(z/λ)
   p₀ + 0.5*(B₀^2 - b^2)
end

ux0(x,z,t) = 0.0
uz0(x,z,t) = 0.0

Bx0(x,z,t) = B₀*tanh(z/λ) + ψ₀*(-π/Lz)*cos(2π*x/Lx)*sin(π*z/Lz)
Bz0(x,z,t) = 0.0 + ψ₀*(-2π/Lx)*sin(2π*x/Lx)*cos(π*z/Lz)

Ay0(x,z,t) = B₀*λ*log(cosh(z))

eq = [
   Dt(ρ(x,z,t)) ~
      -ux(x,z,t)*Dx(ρ(x,z,t)) - uz(x,z,t)*Dz(ρ(x,z,t)) -
      ρ(x,z,t)*(Dx(ux(x,z,t)) + Dz(uz(x,z,t))),
   Dt(p(x,z,t)) ~
      -ux(x,z,t)*Dx(p(x,z,t)) - uz(x,z,t)*Dz(p(x,z,t)) -
      γ*p(x,z,t)*(Dx(ux(x,z,t)) + Dz(uz(x,z,t))),
   Dt(ux(x,z,t)) ~
      -ux(x,z,t)*Dx(ux(x,z,t)) - uz(x,z,t)*Dz(ux(x,z,t)) -
      1/ρ(x,z,t)*(p(x,z,t) + 0.5*(Bx(z,x,t)^2 + Bz(x,z,t)^2)) +
      1/ρ(x,z,t)*(Bx(x,z,t)*Dx(Bx(x,z,t)) + Bz(x,z,t)*Dz(Bx(x,z,t))) +
      ν/ρ(x,z,t)*∇²(ux(x,z,t)),
   Dt(uz(x,z,t)) ~
      -ux(x,z,t)*Dx(uz(x,z,t)) - uz(x,z,t)*Dz(uz(x,z,t)) -
      1/ρ(x,z,t)*(p(x,z,t) + 0.5*(Bx(z,x,t)^2 + Bz(x,z,t)^2)) +
      1/ρ(x,z,t)*(Bx(x,z,t)*Dx(Bz(x,z,t)) + Bz(x,z,t)*Dz(Bz(x,z,t))) +
      ν/ρ(x,z,t)*∇²(uz(x,z,t)),
   Dt(Ay(x,z,t)) ~
      -ux(x,z,t)*Dx(Ay(x,z,t)) - uz(x,z,t)*Dz(Ay(x,z,t)) +
      η*∇²(Ay(x,z,t)),
   Bx(x,z,t) ~ -Dz(Ay(x,z,t)),
   Bz(x,z,t) ~ Dx(Ay(x,z,t))
   ]

domains = [x ∈ Interval(x_min, x_max),
           z ∈ Interval(z_min, z_max),
           t ∈ Interval(t_min, t_max)]

# BCs: periodic in x, Neumann in z
# ICs: set from functions
bcs = [ρ(x,z,0) ~ ρ0(x,z,0),
       ρ(x_min,z,t) ~ ρ(x_max,z,t),
       Dz(ρ(x,z_min,t)) ~ 0.0,
       Dz(ρ(x,z_max,t)) ~ 0.0,

       p(x,z,0) ~ p0(x,z,0),
       p(x_min,z,t) ~ p(x_max,z,t),
       Dz(p(x,z_min,t)) ~ 0.0,
       Dz(p(x,z_max,t)) ~ 0.0,
       
       ux(x,z,0) ~ ux0(x,z,0),
       ux(x_min,z,t) ~ ux(x_max,z,t),
       Dz(ux(x,z_min,t)) ~ 0.0,
       Dz(ux(x,z_max,t)) ~ 0.0,
       
       uz(x,z,0) ~ uz0(x,z,0),
       uz(x_min,z,t) ~ uz(x_max,z,t),
       Dz(uz(x,z_min,t)) ~ 0.0,
       Dz(uz(x,z_max,t)) ~ 0.0,
       
       Ay(x,z,0) ~ Ay0(x,z,0),
       Ay(x_min,z,t) ~ Ay(x_max,z,t),
       Dz(Ay(x,z_min,t)) ~ 0.0,
       Dz(Ay(x,z_max,t)) ~ 0.0,

       Bx(x,z,0) ~ Bx0(x,z,0),
       Bx(x_min,z,t) ~ Bx(x_max,z,t),
       Dz(Bx(x,z_min,t)) ~ 0.0,
       Dz(Bx(x,z_max,t)) ~ 0.0,
       
       Bz(x,z,0) ~ Bz0(x,z,0),
       Bz(x_min,z,t) ~ Bz(x_max,z,t),
       Dz(Bz(x,z_min,t)) ~ 0.0,
       Dz(Bz(x,z_max,t)) ~ 0.0,
      ] 

@named pdesys = PDESystem(eq, bcs, domains,
   [x,z,t],
   [ρ(x,z,t), p(x,z,t), ux(x,z,t), uz(x,z,t), Ay(x,z,t), Bx(x,z,t), Bz(x,z,t)])


## Discretization

order = 2

discretization = MOLFiniteDifference([x=>dx, z=>dz], t, approx_order=order, grid_align=center_align)

## Convert the PDE problem into an ODE problem
println("Discretization:")
@time prob = discretize(pdesys,discretization)

The error messages are

ERROR: LoadError: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 9 and 5")
Stacktrace:
  [1] _bcs1
    @ ./broadcast.jl:516 [inlined]
  [2] _bcs(shape::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, newshape::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}})
    @ Base.Broadcast ./broadcast.jl:510
  [3] broadcast_shape
    @ ./broadcast.jl:504 [inlined]
  [4] combine_axes
    @ ./broadcast.jl:499 [inlined]
  [5] instantiate
    @ ./broadcast.jl:281 [inlined]
  [6] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, Type{Pair}, Tuple{Matrix{Num}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(substitute), Tuple{Tuple{SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}}, Matrix{Vector{Pair{SymbolicUtils.Symbolic{Real}, Num}}}}}}})
    @ Base.Broadcast ./broadcast.jl:860
  [7] BoundaryHandler(bcs::Vector{Equation}, s::MethodOfLines.DiscreteSpace{2, 8, MethodOfLines.CenterAlignedGrid}, depvar_ops::Vector{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}}, tspan::Tuple{Float64, Float64}, derivweights::MethodOfLines.DifferentialDiscretizer{Any, Dict{Function, DiffEqOperators.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/MxvKV/src/bcs/boundary_types.jl:167
  [8] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/MxvKV/src/discretization/MOL_discretization.jl:48
  [9] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/MxvKV/src/discretization/MOL_discretization.jl:146

What does it mean? Which part of the system is not properly set?

[xtalax]

This is the largest system I've seen someone try to solve, which system is it? Let me poke about in the source to see whats causing this

[xtalax]

I notice that you have Bx(z,x,t) as well as Bx(x,z,t), please make these conform to the same argument signature. This change was enough to get it to run on master, I have not tried on v0.3.1 but expect this to be the culprit.

Let me know how you get on!

[henry2004y]

That's a good observation! The discretization step looks ok.

I am at the step of extracting results now. If I use grid size 32*32, this step

solBx = map(d -> sol[d][end], grid[Bx(x, z, t)])

is tediously long. I've waited for more than 30 min, and it's still runnning. I will try a smaller grid size.

[xtalax]

There is also a possibility to use reshape, but you have to calculate grid size:

x = grid[x]
z = grid[z]
Nx = length(x)
Nz = length(z)
@variables Bx[1:Nx,1:Nz](t)
solBx = reshape([sol[Bx[(i-1)*Nz+j]][end] for i in 1:Nx for j in 1:Nz],(Nx,Nz))

Though am unsure if this is faster. This is a serious performance bottleneck that we will be trying to address soon.

[henry2004y]

Thanks again for identifying the mistakes in constructing the equations! I find it error-prone to express a moderately complicated system like this with symbolic codes, and hope the syntax can somehow be improved in the near future. I tried to fix some errors in the equations posted earlier as well as the initial conditions, but I encountered some other problems. I may open up another issue for discussion. This one can be closed if appropriate.

[xtalax]

@henry2004y This bottleneck will be solved in this coming update: https://github.com/SciML/MethodOfLines.jl/pull/154

[xtalax]

@henry2004y just to let you know, the bottleneck should now be closed