MethodOfLines.jl icon indicating copy to clipboard operation
MethodOfLines.jl copied to clipboard

JuliaSimCompiler doesn't work for equation system that works fine with default compiler

Open ctessum opened this issue 1 year ago • 7 comments

Hello,

As suggested in #311, I am experimenting with using the JuliaSimCompiler to allow the use of simulations with larger point counts. I am working with a system similar to this example:

using EarthSciData, EarthSciMLBase, GasChem,
    DomainSets, ModelingToolkit, MethodOfLines,
    Dates, Latexify

@parameters t lev lon lat

start, finish = Dates.datetime2unix.([
    Dates.DateTime(2022, 5, 1),
    Dates.DateTime(2022, 5, 3),
])

domain = DomainInfo(
    partialderivatives_lonlat2xymeters,
    constIC(1.0f0, t ∈ Interval(start, finish)),
    zerogradBC(lon ∈ Interval(-130.0f0, -100.0f0))
)

# Skip unit enforcement for now
ModelingToolkit.check_units(eqs...) = nothing

geos = GEOSFP("0.25x0.3125_NA", t; coord_defaults = Dict(:lat => 34.0, :lev => 1))

model = SuperFast(t) + domain + Advection() + geos

render(latexify(equations(get_mtk(model))))

The equation system looks like this (I cropped the long ones):

lagrida_latex_editor

Basically there is advection driven by external wind fields and and chemical reactions.

When I try to discretize it with the normal compiler, it works fine:

discretization = MOLFiniteDifference([lon => 6], t, approx_order=2)
@time prob = discretize(get_mtk(model), discretization) # great success

However, with the JuliaSimCompiler, it fails:

using JuliaSimCompiler
@time prob = discretize(get_mtk(model), discretization)

JuliaSimCompiler: generate_system
ERROR: ExtraVariablesSystemException: The system is unbalanced. There are 121 highest order derivative variables and 120 equations.
More variables than equations, here are the potential extra variable(s):
... elided ...
 superfast₊HO2[3]
 superfast₊ISOP[3]
 superfast₊NO[3]
 Dt(superfast₊O3[3], 1, false)
 superfast₊O3[4]
 superfast₊NO2[4]
 superfast₊ISOP[4]
 superfast₊HO2[4]
 superfast₊NO[4]
... elided ...
Stacktrace:
  [1] error_reporting(state::IRState, bad_idxs::Vector{Int64}, n_highest_vars::Int64, iseqs::Bool, orig_inputs::Tuple{})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/ZmLra/src/structural_transformation/utils.jl:50
  [2] check_consistency(state::IRState, orig_inputs::Tuple{})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/ZmLra/src/structural_transformation/utils.jl:85
  [3] structural_simplify!(state::IRState, io::Nothing; alias_eliminate::Bool, tearing::Bool, check_consistency::Bool, debug::Bool)
    @ JuliaSimCompiler ~/.julia/packages/JuliaSimCompiler/vX0Sq/src/structural_transformation/dae_transformations.jl:87
  [4] structural_simplify! (repeats 2 times)
    @ ~/.julia/packages/JuliaSimCompiler/vX0Sq/src/structural_transformation/dae_transformations.jl:75 [inlined]
  [5] #structural_simplify#84
    @ ~/.julia/packages/JuliaSimCompiler/vX0Sq/src/structural_transformation/dae_transformations.jl:69 [inlined]
  [6] structural_simplify(::IRSystem)
    @ JuliaSimCompiler ~/.julia/packages/JuliaSimCompiler/vX0Sq/src/structural_transformation/dae_transformations.jl:68
  [7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ MethodOfLinesJuliaSimCompilerExt ~/.julia/packages/MethodOfLines/YO1pd/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl:62
  [8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ MethodOfLinesJuliaSimCompilerExt ~/.julia/packages/MethodOfLines/YO1pd/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl:57
  [9] top-level scope
    @ ./timing.jl:273 [inlined]
 [10] top-level scope
    @ ./Untitled-1:0

In the error message above Dt(superfast₊O3[3], 1, false) sticks out, mostly because it's a derivative and not a variable. Let me know if you have any suggestions of how to fix this!

ctessum avatar Sep 20 '23 16:09 ctessum